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1 ,      Introduction 

This    report  describes    a  numerical   analysis    of    heat    transfer   of    a   typical 
jet  vane   configuration  used   for    thrust   vector   control.      The  work  was    carried 
out  under   contract   Nos.    N62271-85-M-0443   and   N62271-86-M-0206 ,    for    the   Naval 
Postgraduate   School. 

The    tasks    to   be   accomplished   under    the   first   contract   were: 
Task    I:      Formulate    the   conservation   equations    of   momentum  energy   for 
two-dimensional,    supersonic   flow  in  geometries    typical   of    thrust  vector 
control   systems. 

Task    II:      Formulate   boundary   conditions    for    these   equations    appropriate 
to    thrust  vector   control   systems. 

The    tasks    to   be   accomplished   under    the  second   contract   were: 
Task   I:  Continue   and   update    the   formulation   of    thrust  vector   control 

geometries    based   on    the    input   from    the   Naval  Weapons    Center 
(NWC). 
Task    II:  Construct    the   computational  model   for    implementation   in    the 

PHOENICS   code,    of    the    thrust  vector   control   geometries    and 
flow  conditions    provided   by   NWC. 
Task    III:  Run    the   PHOENICS   code    for    the   previously   formulated   models. 

Analyze   and    interpret    the   PHOENICS    results    for   surface 
temperature  "and   heat   flux. 
Thrust  vector   control   components   such  as   jet   vanes    and   jet    tabs    are 
exposed    to   high  speed   hot  gases    at    the   exit   of   a  rocket   nozzle. 

Estimation  of    the   heat    transfer   from   the   hot   exhaust  gases    to    the   vane    is 
major   consideration   in   the   correct   design   of   a  vane,    and   its    ability    to 
survive   during   its   mission. 


The  research  work  was  done  under  the  framework  of  the  tasks.  A  brief 
survey  of  what  has  been  done  according  to  the  task  is  given: 
Task  1  (M-0443):   Heat  transfer  modeling  of  thrust  vector  control  vane 
requires  supersonic  compressible  viscous  flow  analysis. 

In  order  to  meet  the  requirements,  the  conservation  differential 
equations  of  mass  momentum  energy  and  the  two  k-e  turbulent  equations  were 
formulated,  and  additional  algebraic  formulas  for  the  relations  between 
pressure  density  and  the  equation  of  state  for  ideal  gas. 

Task  II  (M-0443):   The  physical  dimensions  of  the  flow  field  grid  were 
chosen  and  the  boundary  conditions  for  the  Navier-Stokes ,  energy  and  the  two 
k-q  turbulent  model  equations  were  given. 

Task  I  (M-0206):   Working  on  the  task,  the  actual  configuration  of  a  jet 
vane  that  is  presently  being  tested  at  NWC  has  been  modeled.   The  geometry 
being  used  is  a  wedge  which  has  the  same  half  angle  and  dimensions  as  the  NWC 
jet  vane. 

Task  II  (M-0206):   BFC  (Body  fitted  coordinate)  version  of  PHOENICS  code 
(Ref.  3)  was  used  for  calculating  the  flow-field  and  heat  transfer  over  the 
model.   Using  the  BFC,  a  better  geometrical  approximation  to  vane  shape  could 
be  achieved. 

Non-Uniform  grids  have  been  utilized  in  order  to  model  complicated  regions 
in  the  flow  field.   Relaxation  parameters  and  false  times teps  options  were 
adjusted  to  enable  efficient  computer  runs  with  good  convergence. 
Task  III  (M-0206):   In  carrying  out  this  task,  four  major  runs  have  been 
analyzed : 

Two  geometric  configurations  were  used:   wedge  vane  and  blunt  vane  (see 
Figures  1,  2,  3,  4);  each  one  in  both  laminar  and  turbulent  flow  conditions. 


Numerical  results  for  fluid  field  and  thermodynamic  properties  of 

pressure,  temperature,  density,  Mach  number  and  velocities  are  given  in 
appendix  C. 

Post-calculations  of  heat  transfer  coefficient,  skin  friction  coefficient 

and  Stanton  number  are  given  in  Figures  (6,  7,  8,  9,  10,  11). 

The  next  chapters  describe  in  more  detail  the  process  of  building  the 

model  and  the  analysis  of  the  results. 


2.   PHOENICS  Description 

The  present  work  addresses  the  heat  transfer  modeling  of  thrust  vector 
control  systems.   In  this  effort  the  Navier-Stokes  approach  is  applied  by 
using  a  computer  code  which  is  capable  of  simulating  a  large  number  of  fluid 
flow,  heat  transfer  and  chemical  reaction  processes  which  arise  in  industry 
and  elsewhere.   This  code  is  called  PHOENICS,  which  is  an  acronym  standing 
for:   'Parabolic,  Hyperbolic  or  Elliptic  Numerical  Integration  Code  Series.1 
The  name  comes  from  the  fact  that  the  differential  equations  of  fluid  flow, 
etc.  arise  in  forms  classified  by  mathematicians  as  parabolic,  hyperbolic  or 
elliptic;  and  PHOENICS  solves  these  equations,  whatever  their  form. 

Built  into  PHOENICS  are  the  major  conservation  laws  of  physics  (mass, 
momentum,  and  energy)  applied  to  a  large  number  of  continuous  subdomains 
called  'cells,'  into  which  the  domain  of  study  is  artificially  divided.   The 
number  of  cells  can  be  few  or  many  according  to  the  requirements  of  the 
problem.   Because  of  numerical  stability  considerations  the  restrictions  on 
cell  refinement  can  become  particularly  burdensome  in  the  calculation  of  a 
turbulent  boundary  layer  where  a  very  fine  mesh  near  the  wall  may  be 
required. 

When  supplied  with  appropriate  information  concerning:   the  physical 
properties  of  the  materials,  the  geometrical  and  other  constraints,  the  inlet 
and/or  initial  conditions,  PHOENICS  computes  the  corresponding  solutions  to 
the  relevant  differential  equations,  expressing  them  as  tables  of  numbers 
describing  the  field  of  velocity,  temperature  concentration,  etc. 

Detailed  information  about  PHOENICS  is  given  in  [Ref.  3]. 


2.1   The  Structural  Principle  of  PHOENICS 

The  code  consists  of  three  major  parts:   Satellite  subroutine, 
Ground  subroutine  and  Earth  library. 

The  satellite  subroutine  is  the  main  input  subroutine  and  should 
provide  the  answers  to  the  questions: 

-  what  kind  of  process  is  to  be  simulated 

-  what  are  the  properties  of  the  fluid 

-  what  are  the  shape  and  size  of  the  domain 

-  how  fine  is  the  grid  to  be  employed 

-  to  what  degree  of  accuracy  is  the  calculation  to  be  continued 

-  and  what  output  should  be  provided 

Ground  subroutine  is  active  during  the  computing  process  and  is  used  for 
updating  properties  which  vary  with  time,  temperature,  etc.   For  example: 
viscosity  depends  on  temperature  or  density  depends  upon  pressure  and 
temperatures,  etc. 

Earth  library  is  the  main  solver  generator.   It  is  given  as  a  binary 
library  and  does  not  enable  the  user  access  to  the  source  code. 

2.2  Numerical  Scheme 

The  numerical  scheme  used  by  the  code  is  the  simpler  (semi-implicit 
method  for  pressure-linked  equations  revised)  (Ref.  9).   The  scheme  was 
developed  by  Patankar,  S.  V.  and  Spalding,  D.  B. 

The  scheme  requires  an  additional  dependent  variable,  the  pressure 
correction,  which  has  no  physical  meaning  but  should  take  part  in  the 
process. 

The  value  of  the  pressure  correction  should  tend  to  zero  in  the 
convergence  process. 

Two  additional  differential  equations  are  solved:   for  the  pressure,  and 
for  the  pressure  correction. 

5 


3.   Geometry  and  Dimensions 

Symmetrical  2-D  planar  geometry,  which  is  shown  in  Figure  2,  was  chosen 
to  be  the  approximation  of  the  MWC  vane  in  Figure  1 . 

Two  geometrical  profiles  were  examined,  one  with  wedge  leading  edge  and 
the  second  with  blunt  leading  edge. 

The  dimensions  of  the  domain  in  Figure  3  and  4  satisfy  aspect  ratio  of 
10:1  in  the  vertical  y  coordinate.   A  high  aspect  ratio  in  the  coordinate  is 
important  for  the  assumption  of  free  stream  conditions  at  the  upper  boundary. 


4.  Assumptions 

Postulating  the  right  or  the  wrong  assumptions  has  the  most  influence  on 
modeling  process.   The  stage  was  carried  out  very  carefully  in  order  to  make 
the  most  compatible  model  with  reality. 

4.1  Steady  state; 

The  modeling  assumes  steady  state  physical  phenomenon  process. 

— —  (all  properties)  =  0 

0  t 

This  is  a  valid  assumption  since  the  time  constant  for  the 
convection  process  is  much  shorter  than  the  time  constant  for  the  wall 
conduction. 

By  assuming  the  wall  temperature  to  be  constant,  the  two  procedures 
are  decoupled. 

In  hot  flow  it  is  important  to  run  the  code  for  a  wide  range  of  wall 
temperature  which  will  take  into  account  the  influence  of  different 
temperatures  on  the  heat  convection  process. 

4.2  Cold  Air  Flow 

Ambient  temperature  air  flow  which  was  utilized  by  NWC  experiments 
is  being  used  in  the  computations. 

4.3  Ideal  Gas 

The  gas  is  assumed  to  satisfy  the  ideal  gas  equation  of  state 
p  =  pRT  (4.1) 

This  is  a  fairly  good  assumption  for  nonreactive  gas  flow.   In  spite  of  the 
values  of  static  temperature  can  decrease  to  200[k],  the  density  remain 
relatively  low. 


This  assumption  is  an  important  simplification  to  the  solution  in  Ref.  10 
which  used  the  isentropic  relation  between  pressure  and  density  instead 

_P_a  (-L.)1^  (4'2) 

Po     po 

4.4  Constant  Pr,^ 

Prantdl  number  and  y  (ratio  of  specific  heats)  were  found  to  have 
negligible  variations  in  the  temperature  range  of  the  model.   (200k  *  350k) 

4 .5  Varying  Viscosity  and  Thermal  Conductivity; 

1j  and  k  are  much  more  dependent  on  temperature  especially  very  close  to 
the  solid  wall  where  values  of  u  and  k  influence  strongly  the  shear  and  heat 
transfer  mechanism.   To  account  for  the  temperature  dependence  power  law 

relations  have  been  formulated  for  y  and  k. 

T   0.666 

M  =  ii0(Y~)  (4-3) 

xo 

X   0.666 
k  =  k  (,**-)  (4.4) 

4.6  Parallel  Flow 

Gas  flow  at  the  exit  of  the  exhaust  nozzle  is  more  likely  to  be  a  conic 

source  flow  than  parallel  flow. 

If  the  half  angle  of  the  nozzle  is  small,  (a  <  15°),  parallel  flow  is  a 
good  assumption 

4.7  Negligible  Radiation 

Assessments  that  were  done  showed  that  heat  convection  is  at  least  one 
order  of  magnitude  greater  than  heat  flux  by  radiation. 


4.8  Laminar  and  Turbulent  Solutions 

In  order  to  overcome  lack  of  ability  to  predict  transition,  separated 
laminar  and  turbulent  calculations  were  done  for  each  case.  The  turbulent 
solution  utilizes  the  (k-e)  eddy  viscosity  model  Ref.  5. 

4.9  Constant  Wall  Temperature 

The  vane  wall  is  assumed  to  have  constant  temperature  during  the  time  of 
calculation. 


5.   Governing  Equations 

The  conservation  equations  for  the  compressible  flow  of  the  mathematical 
model  consists  of  a  viscous,  Newtonian  perfect  gas  consisting  of  the  following 
six  differential  equations: 

Conservation  of  Mass: 

3 1 

Conservation   of   momentum: 

-rr-Cp*)   +  V    (pV>  -    mV<J>)    VP  (5.2) 

where  f  is  V  or  W  velocity  component  for  y  and  z  direction. 
Conservation  of  Energy 

■jI-  (ph)  +  V'Cp^h  -Jj-  Vh)  -J£  (5.3) 

where   h  is    the    total   enthalpy. 

h  =   Cp  TQ 
where   To   is    the    total    temperature 

To  =   Tstat*(l   +  ^M2) 

In  the  case  of  laminar  flow  the  governing  equations  (5.1),  (5.2),  (5.3) 
are  sufficient  to  determine  a  solution  when  proper  boundary  conditions  are 
applied  and  the  equation  of  state  (4.1)  is  provided. 

Turbulence  Model: 

In  turbulent  flow  it  is  necessary  to  hypothesize  a  turbulence   :>del 
relating  the  turbulent  viscosity  to  the  other  problem  variables. 
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The  model  used  in  PHOENICS  is  the  eddy  viscosity  (k-e)  model  [Ref.  3, 
Ref.  5J.   In  this  model  k,  the  turbulent  kinetic  energy  and  e,  the  turbulence 
dissipation  rate,  are  treated  as  properties  of  the  flow  and  conservation 
equations  are  postulated  for  these  properties.   The  two  conservation  equations 
are:   one  for  k,  the  kinetic  energy  of  turbulence: 

£!S  „  _L.  (JS«  *)  +Gk-  e  (5.4) 

Dt   3Xj  Vak   3Xjy    K 

Second  equation  for  e,  the  dissipation  rate  of  turbulence 


De    3  ,veff  3e  >   _e 
Dt  =  dXy    ae  3Xj;   K 


"ST  "  T77T(— —  ^7")  +  -  (clGk    C2e'  (5.5) 


where 

3Ui   3Uj   3UJ 

ueff  "  "lam  +  PSk2/e  (5,7) 

cl»  c2»  ak»  °e»  cu  are  emPiri-cal  constants  which  are  provided  in  PHOENICS. 

The  reason  for  using  the  (k-e)  model  is  because  it  is  the  most  verified 
model  for  engineering  applications.   It  combines  simplicity,  universality,  and 
realism  of  predictions  in  most  cases. 

Two  additional  differential  equations  are  solved  also  in  order  to  satisfy 
the  SIMPLER  algorithm  as  was  mentioned  in  chapter  2.2.  The  description  of  the 
pressure  and  pressure  correction  equations  is  provided  by  Ref.  9. 
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6.   Input  Variables 

The  properties  of  mach  no.  stagnation  presence  and  temperature  of  the  gas 
were  provided  by  NWC;  additional  properties  were  taken  from  air  tables: 

Mach  number: 

M»  "  3.2 

Stagnation  pressure: 

P0  =  55. 105  [Pa] 

Stagnation    temperature: 

To    =    555.55    [K] 

Gas    constant 

R  =   287.    [J/kg'k] 

Specific  heat  ratio 

y   =  1.35 

Laminar  Prandtl  Number 

Pr  =  0.7 

Turbulent  Prandtl  Number 

Prt  =  0.9 

Constant   Pressure   Specific   Heat 

Cp  =   R/(1-1/y)    [J/kg-k] 

Laminar  Viscosity 

u  =  0.1716   *   10"5   *   (T/273.)0-666 

Thermal  Conductivity 

k  -  u  Cp/Pr 

The  gas  properties  in  the  inlet  boundary  are  equivalent  to  the  properties 
at  nozzle  exit.  Inlet  properties  are  calculated  from  the  stagnation  values  in 
the  combustion  chamber.   The  calculation  was  done  by  assuming  one  dimensional 
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isentropic   expansion   from  combustion   chamber    to    the   nozzle   exit    (inlet   for 
the   vane). 


Pressure  P±  =   P0/(l   +  -2~  H2)Y/Y"1 


r-1 
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(6.1) 


Temperature  Tj.  -   TQ/(1   +  JLli  m2)  (6>2) 

Density  Pi  =   P/RT  (6.3) 

Enthalpy  ht  =   CpTi  (6.4) 


Sonic  Velocity        CA  =   /yRTi  (6.5) 

Velocity  Wi  -   OM  (6.6) 

The   subscript    i  signifies    inlet   property 


7 .   Boundary  Conditions 

The  flow  field  described  in  Figures  3,  4  has  four  boundaries,  which  can 
be  named:   inlet,  outlet,  freestream  boundary  and  solid  wall. 

Super  sonic  flows  have  a  hyperbolic  mathematical  nature.   The  field 
consists  of  influence  zones,  the  flow  at  every  point  is  governed  only  by  its 
influence  zone,  basically  by  the  upwind  stream. 

As  a  consequence  from  the  discussion,  it's  obvious  that  the  outlet 
boundary  condition  has  no  influence  on  the  upstream  flow.   The  boundary  values 
that  are  given  at  the  outlet  are  to  satisfy  some  numerical  needs  only. 

7.1  Inlet 

Parallel  uniform  flow  with  known  velocity,  enthalpy,  pressure,  and 
density:   equation  (6.1),  (6.3),  (6.4),  (6.6)  are  given  at  the  left  boundary 
of  the  grid.   In  PHOENICS  this  is  specified  as  the  LOW  side  of  the  first  Z 
cell. 

In  turbulent  flow,  boundary  conditions  are  supplied  for  k  and  e.   The 
values  that  are  given  are  based  on  empirical  values: 

kt  =  0.0  V^  (7.1) 

ei  =   0.16   kl.5/(5*GH)  (7.2) 

where   GH   is    half    the   vane    thickness 

7.2  Outlet 

As  was  mentioned  previously,  the  outlet  has  negligible  effect  on  the 
results.   The  only  property  that  is  specified  at  the  outlet  is  the  pressure. 
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7.3  Frees tream  Boundary 

Assuming  that  the  upper  boundary  is  chosen  to  be  far  enough  away,  the 
default  boundary  condition  option  of  PHOENICS  is  used.   This  implies  a  line  of 
symmetry  where  all  gradients  are  zero. 

7.4  Solid  Wall 

Zero  velocity  and  constant  wall  enthalpy  (temperature)  are  assumed  on  the 
wall.   In  PHOENICS  the  wall  is  the  SOUTH  side  of  the  first  y  cell.   The  high 
enthalpy  and  velocity  gradients  near  the  wall  demands  a  refined  grid  close  to 
the  wall.   Values  of  shear  stress  and  heat  flux  are  calculated  to  first  order 
accuracy  using: 

3W       Wl  n    n 

tw  =  u  —  =  u  - — t-7?  t  / .  l ; 

w    8y    Ayi/2 

qw   Pr  3y  '  Pr   Ayi/2  U'/; 

In  turbulent  flow,  a  wall  function  is  used  to  provide  the  wall  condition 
for  velocity,  enthalpy,  k,  and  e 

7.5  Wall  Function 

The  wall  problem  in  the  numerical  computation  of  flows,  especially  in 
turbulent  flow,  is  an  old  one  and  most  authors  have  adopted  similar 
techniques.   In  effect  they  "bridge  over"  the  region  very  close  to  the  wall  by 
introducing  special  functions  which  are  called  wall  functions.   These  are 
often  empirical  in  origin.  Accounts  may  be  found  in  Ref.  11. 

The  problem  arises  as  follows.   Turbulence  dies  out,  close  to  the  wall, 
because  the  no  slip  condition  and  the  rigidity  of  the  wall  make  all  the 
velocity  components  fall  to  zero.   The  consequence  is  that  the  effective 
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viscosity  and  other  transport  properties  fall  there  to  their  laminar  values 
and  the  result  is  a  rapid  variation  with  distance  from  the  wall  both  of  the 
<j>  's  and  of  their  gradients. 

Where  <f>  signifies  general  dependent  variable,  it  is  possible  to  compute 
these  variations  in  detail,  by  using  a  computer  code  such  as  PHOENICS  on  two 
conditions: 

(i)   the  grid  points  must  be  packed  into  the  region  of  steep  gradient 
changes  closely  enough  for  sufficient  numerical  accuracy  to  be  obtained 

(ii)   the  functions  appearing  in  the  turbulence  model  equations  must 
properly  represent  the  influence  of  local  Reynolds  number  on  turbulence. 

Under  the  conditions  above,  the  wall  function  sequences  in  the  program 
act  as  follows: 

The  Reynolds  number  is  first  evaluated,  based  on  the  resultant  velocity 
parallel  to  the  wall,  on  the  distance  from  the  wall  to  the  grid  node  and  on 
density  and  laminar  viscosity.   If  this  Reynolds  number  is  less  than  132.25 
(the  value  at  which  the  laminar  and  turbulent  wall  function  intersect)  a 
laminar  wall  function  is  used.   If  this  Reynolds  number  turns  out  to  be 
greater  than  132.25  the  velocity  variation  is  logarithmic  and  the 
corresponding  shear  stress  coefficient  is  evaluated.   This  corresponds  to  the 
commonly  used  "log  law"  wall  function.   [Ref.  4] 

t 

7 .6      Boundary  Conditions    in   Phoenics 

PHOENICS   utilizes    source    terms    for   creating   boundary   conditions.      The 
form  of    the   source    term  of   each  dependent   variable   <J>   is: 

S*  =   ([m]    +  C(())    CVcf.  -   <j,p)  (7.3) 

where:      m  -   is   mass    flux  source 

4>p  -   is    the   value   of    the   dependent   variable   at   point  near    the   boundary 

C<J),    V<J>  -   two   coefficients   specified   by    the   user.      The   source    term   for 
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mass  flux  is  simply 

Sm  "   Cm   (Vm  -   Pp)  (7.4) 

where:      Pp  -is    the   pressure   near    the   boundary   and   C^,    Vm  are    two   coefficients. 
The   values    of   C<J>  and   V<j>  for    the   dependent  variables    in   SATELLITE   are:      At    the 
Inlet: 

c*=2FrW  (7-5) 

vm  =   po  Pi/Po  (7.6) 

Cw  -   Ch  =   Ck  =   C£  =   0.  (7.7) 

(7.8) 
(7.9) 
(7.10) 
(7.11) 

(7.12) 

Vm  =   Pi  (7.13) 

At    the  Wall    (laminar) 

Cw  =   u/(0.5Au!)  (7.14) 

Vw  =  °  (7.15) 

Ch  =   u/Pr/(0.5-Ayi)  (7.16) 

Vh  -   S*Tw  (7-17) 

At    the  Wall    (turbulent) 

cw  "   ch  =   ck  =   ce  =  WALL  (7.18) 

Vw  -  Vk  -  V£  =  0  (7.19) 

Vh  »   Cp*Tw  (7.20) 


Vw  = 

wi 

Vh  = 

hi 

Vk   = 

Ki 

v£  = 

ei 

At 

the   Outlet: 

Cm  = 

1000 

* 

wi 

Wpi 
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8.  Mesh  Generation 

In  this  work  a  two-dimensional  mesh  is  being  used  with  18  x  29  cells  in 
the  y  and  z  coordinate  respectively.   A  Nonuniform  grid  has  been  used  for 
both  directions.   Figures  3  and  4  shows  the  grid  in  the  z  direction.   A  finer 
grid  is  used  in  the  blunt  region,  IZ  ■  (7  4  17),  and  in  the  zone,  where  the 
inclined  wall  transitions  to  a  straight  wall,  IZ  =  (23  *  26). 

In  the  y  coordinate,  except  in  the  boundary  layer  region,  the  grid  is 
uniform.   To  obtain  a  finer  grid  resolution  in  the  boundary  layer  for  the 
laminar  flow  case  the  first  five  cells  in  the  y  direction  from  the  wall  obey 
the  following  proportionality  relationship: 

BYFRAC  (IY)  =  (-f^)3  (f^f1)  (8.1) 

Where  BYFRAC( IY)  is  the  distance  from  the  south  side  to  the  north  side  of  the 
cell  of  particular  interest,  divided  by  total  length  of  the  domain,  IY  is  the 
cell  number,  Amax  is  maximum  allowable  cell  height,  and  GH  is  the  half 
thickness  of  the  TVC  jet  vane. 

A  fine  grid  resolution  for  the  turbulent  flow  case  is  set  up  in  the  same 
way  as  laminar  flow.   The  only  difference  comes  from  the  selection  of  the 
first  five  cells  in  y  direction.   The  following  calculation  shows  the 
difference. 

From  the  laminar  solution  and  the  given  properties  the  following  are 
known: 

w  »  885.2[m/s] 

Ulam  =  l*10"5  (N.s/m] 

Po  =  5.5*106[Pa] 

Pstatic   =   1.048*105[Pa] 

Y  =    1.35 
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p  =  1.835  Ikg/m] 

Using  the  values  above  and  the  length  of  vane,  which  is  0.095m,  A 
corresponding  Reynolds  number  was  calculated: 

P«°  Woo  z   (1.835*888.5*0.095)    ,  rA  *  1n6 

Re* =  ^7^     r7^ 1,s      10 

Using   a  power   law  correlation   for    the   boundary   layer    thickness: 

-  =   0.37    *  Rez"l/5  (8-2) 

z 

From  equation  (8.2)  the  boundary  layer  thickness  at  the  high  end  of  the 

domain  has  been  calculated  as  6  «  2  *  10"*3[m] 

Ay 
With  Re  based  on  W„  the  velocity  parallel  to  the  wall,  -?-   the  distance 

2 

from   the  wall    to    the   first  grid    node,    pM   the   density,    and   ulam   the   laminar 

viscosity,    Ay  must  satisfy    the   condition 

m  Poo     co     y   >   132#25  or  AY  >   6.48   *   10"6[m] 

^A  2    yiam 

Therefore    the    interval   of   Ay   is    chosen   such   that 

2    *    10~3[m]    >   AY   >  6.48   *   10"6[m] 
In    this    effort  using    the   relationship 

BYFRAC(IY)    -   (^-)2    (^g§) 

Ay  has  been  calculated  as  Ay  =  8  *  10~5[m]  which  is  in  the  required  interval. 

For  both  the  laminar  and  turbulent  cases,  cells  in  the  z  direction  were 
adjusted  so  that  the  points  where  possible  physical  phenomena  such  as  shock 
waves  and  expansion  fans  are  expected,  very  fine  cells  were  used.   In  the 
other  parts  of  the  domain  larger  cells  were  used. 
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9.   Heat  Transfer  Analysis 

Skin  friction  and  heat  transfer  quantities  were  calculated  in  both 
laminar  and  turbulent  cases  and  they  are  shown  in  Figures  (6  *  11). 

9 .1   Laminar  Calculation 

In  laminar  flow  fluxes  can  be  derived  directly  from  the  gradients  near 
the  wall.   The  first  cell  is  close  "enough"  to  the  wall  and  gradients  of 
velocity  and  enthalpy  do  not  change  much  in  this  region  near  the  wall.   The 
shear  stress  and  heat  flux  in  the  laminar  case  will  be: 
wl 

hi  -  hw 
*"  ~   ?7  AYi/2  (7*2) 

The  skin  friction  coefficient  and  Stanton  number  will  be: 

00  00 

st  =  qw/ [ p-y.( hr  -  M]  (9-2) 

where  hr  is  the  recovery  enthalpy 


h        r(^-}  M2 

"o =  i  +  5ED  7 

2      oo 


(9.3) 


r-   is    the   recovery   factor 

r  »  /Pr     (laminar   flow)  (9.4) 
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The  coefficient  of  heat  transfer  in  convection  was  calculated  using 

hc  -  P.  "co  Cp  St  (9.5) 

9.2  Turbulent  Calculations 

In  turbulent  flow  the  gradients  of  velocity  and  enthalpy  near  the  wall 
are  very  steep  and  change  rapidly  with  distance  from  the  wall. 

Direct  calculation  of  flux  gradients  is  not  accurate  in  this  case.   The 

log  law  approach  is  used  to  calculate  skin  friction.   In  the  calculations 

using  PHOENICS  flow  field,  the  following  relation  has  been  used. 

2  pw  ^w  ,a   ,v 

Cf  =  — (9.6) 

w2Po»3.33 

CD 

To  obtain  equation  9.6,  the  turbulent  kinetic  energy  equation  has  been  used  as 
a  starting  point.   [Ref .  5] , 

,  .  Mt  ,3u2    r      p2  k  . 

The   source    terra  of    the    turbulent   kinetic   energy   equation  should   be   zero   near 
the  wall   which  means 

£(iH)2     -  C^     =   0  (9.8) 

k  3y     w  Vt 

therefore  the  shear  stress  on  the  wall  can  be  defined  as: 

tw  =  CD  I/2  p,,^       ^  (9.9) 

where  l^  is  the  turbulent  kinetic  energy  on  the  wall,  pw  is  the  density  on  the 

wall  and  Cq  =  0.09  [Ref.  5] ,  substituting  the  values  above  into  the  Blasius 

skin  friction  relation  the  Cf  equation  becomes: 
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2   tw         Pw       2  kw 

cf  = ~  3  « r  "TTIT 


w 


2         Pa 


(9.10) 


w 


The   heat    transfer   quantities    are   evaluated   from    the   Chilton-Colburn   form  of 
Reynolds    analogy. 

st  =   (Cf/2)    *   Pr"2/3  (9.11) 

qw  -  st*P-*u-*<hr-V  (9.12) 

where  equation  (9.3)  is  used  to  evaluate  hr  with  the  recovery  factor  given  as: 

1/3 

r  =  Pr   (turbulent  flow)  (9.13) 

The  convective  heat  transfer  coefficient  is  calculated  by  using  equation  (9.5) 
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10.   Code  and  Computer 

PHOENICS  81,  Body  Fitted  Coordinate  (BFC)  version  has  been  used  in  the 
computations  (see  Ref.  3).   PHOENICS  has  been  installed  on  NPS  IBM  3033  MVS 
1.3  computer.   400  sweeps  per  computer  run  provided  a  reasonable  convergence 
in  all  runs  except  the  turbulent  blunt  case  continuity  error  of  less  than 
4.10~t+  has  been  achieved  in  the  three  runs. 

The  continuity  error  is  the  total  summation  of  the  absolute  mass 
imbalance  in  all  cells  divided  by  the  inlet  mass  flux.   CPU  time  consumption 
varies  from  case  to  case  as  follows: 

Laminar  Wedge      630    CPU  Seconds 

Turbulent  Wedge    630    CPU  Seconds 

Laminar  Blunt      630    CPU  Seconds 

Turbulent  Blunt    1542   CPU  Seconds  for  1000  sweeps 
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11 .   Results  and  Discussion 

The  results  of  the  calculations  are  available  on  appendix  c.   The  tabular 
results  include  the  values  of  pressure,  velocities,  enthalpy,  temperature  mach 
number,  density,  turbulent  kinetic  energy  and  rate  of  turbulent  dissipation. 
The  values  are  given  in  18  x  29  cells  points. 

Skin  friction  and  heat  transfer  results  are  shown  in  Figures  (5-11). 
Laminar  and  turbulent  skin  friction  and  Stanton  number  in  wedge  flow  show 
improvement  compared  to  the  results  reported  by  Yukselen  (Ref.  10).   The  lines 
are  smoother  and  the  oscillations  at  the  end  were  eliminated.   Basically  the 
magnitudes  are  similar  to  those  in  Ref.  10. 

Laminar  blunt  values  are  similar  except  near  the  beginning.   The 
beginning,  as  expected  in  blunt  zone,  creates  higher  rates  of  heat  transfer. 
Even  though  the  blunt  geometry  used  is  a  multi-wedge  shape  it  should  predict 
the  correct  values  except  for  the  stagnation  point  itself. 

Turbulent  blunt  skin  friction  has  different  behavior.   It  has  a  very 
large  value  at  the  first  point  and  then  undershoots  to  values  that  are  smaller 
than  for  wedge.   It  should  also  be  kept  in  mind  that  the  convergence  of  this 
case  wasn't  very  successful. 
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12.   Conclusions  and  Recommendations 

1.  PHOENICS  was  found  to  be  a  friendly  code  for  simulating  complicated 
mixed  heat  transfer  fluid  dynamics  problems. 

2.  Derivation  of  heat  transfer  properties  to  a  vane  solid  wall  in 
laminar  and  turbulent  flow  has  been  installed  in  the  code.   It  can  be  used  for 
predictions  of  heat  transfer  rate  in  both  cold  and  hot  gas  flow. 

3.  Two  features  have  been  added  to  the  code  in  NPS:   The  restart  option 
and  the  use  of  initial  field,  make  it  possible  to  simulate  time  dependent 
processes  and  solve  the  temperature  variation  in  the  vane  itself. 


25 


LIST  OF  REFERENCES 

1.  Baldwin  and  McCormak,  "Numerical  Solution  of  the  Interaction  of  a  Strong 
Shock  Wave  with  a  Turbulent  Boundary  Layer,"  AIAA  Paper  74-558,  AIAA  7th 
Fluid  and  Plasma  Dynamics  Conference  Palo  Alto,  Calif.,  June  17-19,  1974. 

2.  Shang,  J.  S.,  Hankey,  W.  L. ,  and  Law,  C.  H. ,  "Numerical  Simulation  of 
Shock  Wave-Turbulent  Boundary  Layer  Interaction,"  AIAA  Paper,  76-95,  AIAA 
14th  Aerospace  Sciences  Meeting,  Washington,  D.C.,  Janaury  1976. 

3.  Gunton,  M.  C. ,  Rosten,  H.  L. ,  and  Spalding,  D.  B.,  Phoenics  Instruction 
Manual,  Spring  1983,  CHAM  Co,  London. 

4.  White,  Frank  M. ,  Viscous  Fluid  Flow,  Mc  Graw-Hill,  Inc.  1974. 

5.  Launder,  B. ,  and  Spalding,  D.  B. ,  Lectures  in  Mathematical  Models  of 
Turbulence,  Department  of  Mechanical  Engineering  Imperial  College  of 
Science  and  Technology,  London,  England,  Academic  Press,  1972. 

6.  Shapiro,  Ascher  H. ,  The  Dynamics  and  Thermodynamics  of  Compressible  Fluid 
Flow,  Volume,  Department  of  Mechanical  Engineering,  Massachussetts 
Institute  of  Technology,  Newyork,  The  Rolald  Press  Company. 

7.  Schlichting,  Hermann,  Engineering  University  of  Braunschweig  Newyork  St. 
Louis,  Boundary  Layer  Theory  Mc  Graw  Hill  Book  Company,  San  Francisco 
Toronto,  1968. 

8.  Lin,  C.  C,  Turbulent  Flows  and  Heat  Transfer,  Princeton,  New  Jersey 
Princeton,  New  Jersey  Princeton  University  Press,  1968. 

9.  Patankar,  S.  V.,  Numerical  Heat  Transfer  and  Fluid  Flow,  Mc  Graw  Hill, 
Inc.,  1980. 

10.  Yukselen,  A.,  Heat  Transfer  Modeling  of  Thrust  Vector  Control,  Msc. 
Thesis,  Naval  Postgraduate  School  1986. 

11.  Spalding  D.  B.,  "A  General  Computer  Code  for  Two-Dimensional  Elliptic 
Flows,"  Imperial  College,  London,  1977. 


26 


c 
o 


3 


o 

c 
o 

> 


o 


All  Dimensions  in  Millimeters 


59.0 


■90.8- 


4 30 » 


12.7 


Figure   2:     NWC    Jet    Vane    Approximation 
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Figure  3:  Wedge  vane  domain  and  grid 
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Figure  4:   Blunt  vone  domain  and   grid 
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Appendix  A 
Satellite  Listing 
Two  subroutines  SATELLITE  and  GROUND  had  to  be  changed  and  improved.   The 
full  list  is  enclosed  in  Appendix  A  and  B. 

VAN4SAT  and  VANTSAT  are  the  laminar  and  turbulent  SATELLITE  subroutines, 
the  first  has  the  blunt  geometry  and  the  second  has  the  wedge  (it  can  be 
changed  easily  from  wedge  to  blunt  and  vice  versa)  VAN4GRD  and  VANTGRD  are 
exactly  the  same.   They  are  the  GROUND  subroutines,  VANTGRD  is  given  in 
Appendix  B. 
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FILE:  VANTSAT   FORTRAN   Al 


CSDIRECTIVEXXSATLIT   AMI  LEITNER 

C  LAMINAR  SOLUTION  FOR  NWC5    NY=18  NZ=29  YN=GTH 

C  LECSAT  CONVERTED  TO  DIAMSAT 

C  *FILE  NAME:  MODBFCST.FTN 

C  ^ABSTRACT:  SATELLITE  MODEL  MAIN  PROGRAM.  THIS  VERSION  IS 

C  FOR  USE  WITH  THE  BODY-FITTED  COORDINATE  SCHEME  (SUMMER  1984 

C  VERSION)  PROVIDED  AS  AN  ATTACHMENT  TO  SPRING  1983  PHOENICS. 

C  XDOCUMENTATION:  PHOENICS  INSTRUCTION  MANUAL  (SPRING  1983) 

C  WITH  BODY-FITTED  COORDINATES  INSTRUCTION  SUPPLEMENT 

C  (SUMMER  1984). 

C  ^AUXILIARY  SUBROUTINES  (TAPES,  ETC.)  ARE  IN  SATELLITE  LIBRARY 

C  SERVICEU,  WHICH  MUST  BE  INCLUDED  IN  LINK  EDIT  TO  RUN. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  1  STARTS: 


CHAPTER  1 
C— 


COMMON  BLOCKS  AND  USER'S  DATA 


INCLUDE  (CMNGUS) 
INCLUDE  (CMNGRF) 

INCLUDE  (GUSSEQ) 

COMMON/CPI/IPWRIT, IDUM(243) 

DIMENSION  GDTAPE(3),DFAULT(4) 

DIMENSION  ARRAY1(309),ARRAY2(194),ARRAY3(421) 

LOGICAL  ARRAY1,LSPDA,WRT,RD,NAMLST 

INTEGER  ARRAY2,XPLANE,YPLANE,ZPLANE 

INTEGER  P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS,EP,H1,H2,H3,C1,C2, 
SC3,C4 

REAL  NORTH, LOW 

LOGICAL  BFC 

EQUIVALENCE  ( ARRAY1 ( 1 ), CARTES) , (ARRAY2( 1 ), NX) 

EQUIVALENCE  (ARRAY3( 1 ) , SPAREK 1 ) ) , (Ml , Rl ) , (M2, R2) 

EQUIVALENCE  ( LSTRUN, INTGR( 12) ) , (NAMLST, L0GIC(88) ) 

EQUIVALENCE  (LOGIC(20) , BFC) 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  1  ENDS. 
C$DIRECTIVEX*CMNBF1$$ 

C   THIS  FILE  CONTAINS  SATELLITE  COMMON  BLOCKS  FOR  BFC'S 
C   Fl  MUST  BE  DIMENSIONED  TO  GREATER  THAN  OR  EQUAL  TO 
C  (NX+NY+17XNZ+24XNXXNY+6*(NX+1)X(NY+1)+6XND).  THE  VALUE 
C   OF  THE  DIMENSION  MUST  BE  SET  AS  NBFC  IN  GROUP  6  OF  SATLIT. 

COMMON/FOB/FK5000) 

COMMON/CIB/ND/CIC/KOORD 

COMMON/CID/KDBGG, KDBGMF, KDBGCD, KDBIND, KDBMFX, KDBCDT, KDBPCS, 
8  KDBGUV,KDBGPV 

COMMON/CIE/KDBGS,KDBINS 

COMMON/CIF/IGEN/CIG/NCART 
C   THE  FOLLOWING  ARRAYS  MUST  BE  EXACTLY  DIMENSIONED  FOR  NXP1, 
C   NYP1  AND  NZP1,  BUT  MAY  BE  OVER  DIMENSIONED  FOR  ND. 
C   THE  BFRAC  ARRAYS  MUST  BE  DIMENSIONED  TO  ALLOW  FOR  SETTINGS 
C   IN  SATLIT,  THEY  MAY  BE  OVER  DIMENSIONED. 

C0MM0N/CRA/XW(19,30,1)/CRB/XE(19,30,1) 
8       /CRC/YS(2,30,1)/CRD/YN(2,30,1) 
8       /CRE/ZL(2,19,1)/CRF/ZH(2,19,1) 

8       /CRG/RC0N/CRH/DARCY/CRI/BXFRAC(99)/CRJ/BYFRAC(99) 
8       /CRK/BZFRAC(99) 

C0MM0N/CLA/ST0RSA(6 ) , ST0RWD(6 ) , STORP, STORPE, STORPN, 
8  ST0RPH,ST0R1,ST0R2,ST0R3,ST0UNV,PRTBFC,ST0CRN 

COMMON/CLC/BFPLOT 

LOGICAL  STORP, STORPE, STORPN, STORPH, STOR1 , ST0R2, ST0R3, 
8         STORSA, STORWD, STOUNV, PRTBFC, BFPLOT, STOCRN 
C      END 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  1  STARTS: 
C      GRAFFIC  ARRAYS  DIMENSIONED  AS  NEEDED... 

COMMON/GRAF1/PHI1 ( 1 )  /GRAF2/PHI2( 1 ) 
C    .  POROSITY  8  SPECIAL  DATA  ARRAYS  DIMENSIONED  AS  NEEDED... 

DIMENSION  PE(1,1,1),PN(1,1,1),PH(1,1,1),PC( 1,1,1) 

DIMENSION  LSPDA(1),ISPDA(1),RSPDA(1) 
C      USER  PLACES  HIS  VARIABLES,  ARRAYS,  EQUIVALENCES  ETC.  HERE. 

EQUIVALENCE(RAIR,RE(21)),(GAMA,RE(22)),(GSWP,RE(23)) 
1,(GPR,RE(24)),(TW,RE(25)),(GEMU1,RE(26)),(JEMU1,INTGR(D) 
C      USER  PLACES  HIS  DATA  STATEMENTS  HERE. 

DATA  NLSP,NISP,NRSP/1,1,1/ 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  1  ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  2  STARTS: 


VANOOOIO 
VAN00020 
VAN00030 
VAN00040 
VAN00050 
VAN00060 
VAN00070 
VAN00080 
VAN00090 
VAN00100 
VAN00110 
VAN00120 
VAN00130 
VAN00140 
VAN00150 
VAN0016O 
VAN00170 
VAN00180 
VAN00190 
VAN00200 
VAN00210 
VAN00220 
VAN00230 
VAN00240 
VAN00250 
VAN00260 
VAN00270 
VAN00280 
VAN00290 
VAN00300 
VAN00310 
VAN00320 
VAN00330 
VAN00340 
VAN00350 
VAN00360 
VAN00370 
VAN00380 
VAN00390 
VAN00400 
VAN00410 
VAN00420 
VAN00430 
VAN00440 
VAN00450 
VAN00460 
VAN00470 
VAN00480 
VAN00490 
VAN00500 
VAN00510 
VAN00520 
VAN00530 
VAN00540 
VAN00550 
VAN00560 
VAN00570 
VAN00580 
VAN00590 
VAN00600 
VAN00610 
VAN00620 
VAN00630 
VAN00640 
VAN00650 
VAN00660 
VAN00670 
VAN00680 
VAN00690 
VAN00700 
VAN00710 
VAN00720 
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CHAPTER  2    SET  CONSTANTS,  AND  ARRANGE  FILE  MANIPULATIONS 


C 
C 


CD 
CD 


PLEASE  DO  NOT  ALTER,  OR  RE-SET,  ANY  OF  THE  REMAINING 
STATEMENTS  OF  THIS  CHAPTER. 

DATA  CELL, EAST, WEST, NORTH, SOUTH, HIGH, LOW, VOLUME/ 
8  0.,1.,2.,3.,4.,5.,6.,7.  / 

DATA  P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS,KE,EP,H1,H2,H3,C1,C2, 
SC3,C4/1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20/ 
DATA  FIXFLU, FIXVAL , ONLYMS, WALL/1 . E-10 , 1 . E10, 0 . 0 , -10 . 0/ 
DATA  IPLANE,XPLANE,YPLANE,ZPLANE/0,1,2,3/ 

DATA  WRT,  RD,  DFAULT/  .TRUE. ,  .  FALSE. ,  4HDEFA,  4HULT..,  4HDTA/,  1HG/ 
DATA  GDTAPE/4HGUSI,4HE1.D,2HTA/ 
DATA  NL DATA, NI DATA, NRDATA/309, 194,421/ 
DATA  NLCREG,NTCVRG/60,350/ 

DATA  TITPP,TITC1,TITC2,TITC3/3HRH0,4HMACH,4HTEMP,4HCFST/ 
CALL  TAPES(10,GDTAPE,3,1,4XNRDATA) 

READ  DEFAULT  FILE  IF  BLOCKDATA  ABSENT 

IF(INTGR1(29).NE.10)  GO  TO  2 
CALL  WRIT40(40HDATA  ESTABLISHED  IN  BLOCK  DATA.  ) 

GO  TO  3 
CALL  DEFLT 

TAPES(1,DFAULT,4,2,4*NRDATA) 

DATAIO(RD,l) 

WRIT40(40HDATA  TAKEN  FROM  DEFAULT. DTA  ON  GROUP  A/C) 

WRIT40C40HFILE  MODSTL.FTN  IS  THE  SATLIT  USED.      ) 


CALL 
CALL 
CALL 
CALL 


L0GIC(89)  =  .TR'JE, 


DEFINE  DATA  FOR  NRUN  RUNS. 


CHAPTER  3 

c 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

C—  GROUP   41MULTI-RUNS  i  RUN(  1-30X  .T  . , 
C 

RUN(1)=. FALSE. 
C   NOTE.  ALL  RUNS  ARE  DEACTIVATED  AT  THIS 
C   ====   SWITCH  ON  ONE  ONLY  OF  RUNS  1-4  I 

RUN(1)=.TRUE. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
DO  10  IRUN=1,30 

IF( .NOT.RUN(IRUN))  GO  TO  10 
NRUN=NRUN+1 
LSTRUN=IRUN 
10  CONTINUE 

DO  999  IRUN=1,LSTRUN 

IFC.NOT.RUN(IRUN))  GO  TO  999 
INTGR(ll)  =  IRUN 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C—   ALL  INTEGER  VARIABLES  ARE  DEFAULTED 
C      TO  0.0,  UNLESS  OTHERWISE  INDICATED. 
C      E.G.  BY  VARIABLE<10>,  OR  <10.0>  AS 
C     THE  DEFAULT  SETTINGS  OF  ALL  LOGICAL 
C      INDICATED,  E.G.  VARIABLE< .T . >,  OR 
C 

C RUN1 

c 

C GROUP  1.  FLOW  TYPE  i 

C     PARAB< . F . > , CARTES< . T . > , ONEPHS< . T . > 
c 

C—   GROUP  2.  TRANSIENCE  . 


STANDARD  SECTION  2  ENDS, 
USER  SECTION  2  STARTS: 
29*.F.> 


POINT  -  USER  SHOULD 
N  NEXT  STATEMENT. 

USER  SECTION  2  ENDS. 
STANDARD  SECTION  3  STARTS. 


STANDARD  SECTION  3  ENDS. 
USER  SECTION  3  STARTS: 
TO  0,  AND  REAL  VARIABLES 

APPROPRIATE. 

VARIABLES  ARE  ALWAYS 
VARIABLE<.F.>. 


C 

c 
c 
c 

c 

c— 
c 
c 
c 

c 


STEADY<.T.>,ATIME,LSTEP<1>,FSTEP<1> 

TLAST<1.E10>,TFRAC(1-30X30X1.> 

SERVICE  SUBROUTINE  FOR  'NT*  POWER-LAW  TIME  STEPS. 

CALL  GRDPWRCO, NT, TLAST, POWER) 


GROUP  3.  X-DIRECTION  . 
NX<l>,XULAST<1.0>,XFRAC(l-30) 
SERVICE  SUBROUTINE  FOR  POWER-LAW  GR 
CALL  GRDPWRU,  NX,  XULAST,  POWER) 


ID. 


VAN00730 
VAN00740 
VAN00750 
VAN00760 
VAN00770 
VAN00780 
VAN00790 
VAN00800 
VAN00810 
VAN00820 
VAN00830 
VAN00840 
VAN00850 
VAN00860 
VAN00870 
VAN0088O 
VAN00890 
VAN00900 
VAN00910 
VAN00920 
VAN00930 
VAN00940 
VAN00950 
VAN00960 
VAN00970 
VAN00980 
VAN00990 
VAN01000 
VAN01010 
VAN01020 
VAN01030 
VAN01040 
VAN01050 
VAN01060 
VAN01070 
VAN01080 
VAN01090 
VAN01100 
VAN01110 
VAN01120 
VAN01130 
VAN01140 
VAN01150 
VAN01160 
VAN01170 
VAN01180 
VAN01190 
VAN01200 
VAN01210 
VAN01220 
VAN01230 
VAN01240 
VAN01250 
VAN01260 
VAN01270 
VAN01280 
VAN01290 
VAN01300 
VAN01310 
VAN01320 
VAN01330 
VAN01340 
VAN01350 
VAN01360 
VAN01370 
VAN01380 
VAN01390 
VAN01400 
VAN01410 
VAN01420 
VAN01430 
VAN01440 
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C- 

C 

C 

C 

C- 

C- 

C 

C 

C 

C- 

C- 

C 

C 

C- 

C 

C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


GROUP  4.  Y-DIRECTION  < 

NY<1>,YVLAST<1.0>,YFRAC(1-30),RINNER,SNALFA 
SERVICE  SUBROUTINE  FOR  POWER-LAW  GRID* 
CALL  GRDPWRC2, NY, YVLAST, POWER) 
NY=18 


GROUP  5.  Z-DIRECTION  i 

NZ<l>,ZWLAST<1.0>,ZFRAC(l-30) 

SERVICE   SUBROUTINE   FOR   POWER-LAW   GRID: 

CALL    GRDPWR(3,NZ,ZWLAST, POWER) 

NZ=29 


GROUP   6.    MOVING  GRID   OR   DISTORTED    (BODY-FITTED)    GRID    « 

—  MOVING  GRID    : 
MGRID,IZW1,IZW2,AZW2,BZW2,CZW2,PINT,ZW2M1T 

—  BODY-FITTED   GRID   — 
BFC<.T.>,IGEN<1>,ND<1>,NBFC<5000>,KOORD,RCON 
BXFRAC(1-NXK1.0,NXM1X0.0> 
BYFRACC1-NYK1 .  0,  NYM1*0  .  0> 

BZFRACC  1-NZX1 .  0,NZM1*0  .  0> 

SERVICE   SUBROUTINE   FOR   SUB-DOMAIN    SPECIFICATION    (FOR    IGEN=1 

ONLY): 

CALL  DOMAIN( ID, IXF, IXL , IYF, IYL , IZF, IZL ) 
XE(1-NYP1,1-NZP1,1-NDK(NYP1*NZP1XND)*1.0>, 
XW(1-NYP1,1-NZP1,1-ND), 

YN(1-NXP1,1-NZP1,1-NDX(NXP1XNZP1XND)X1.0>, 
YS(1-NXP1,1-NZP1,1-ND), 

ZH(1-NXP1,1-NYP1,1-NDX(NXP1XNYP1XND)X1.0>, 

ZL(1-NXP1,1-NYP1,1-ND),5T0RSA(1-6K6X.F.>,ST0RWD(1-6X6X.F.>, 
STORP< . F . >, STORPE< . F . > , STORPN< . F . > , STORPH< . F . >, STOUNV< . F . >, 
PRTBFC< . F . >, DARCY, BFPLOT< . F . > 

CYCLIC  BOUNDARY  CONDITIONS  ARE  DEFAULTED  INACTIVE  ; 
TO  ACTIVATE  THEM  AT  SELECTED  IZ  SLABS  USE  SERVICE  SUBROUTINE: 

CALL  XCYIZdZ,  .TRUE.) 
SERVICE  SUBROUTINE  TO  DEACTIVATE  CURVATURE  TERMS  IN  U,  V 
AND  W  EQUATIONS  ASSOCIATED  WITH  CURVATURE  OF  IX,  IY,  IZ 
GRID  LINES  RESPECTIVELY: 

CALL  UCURVEdZ,  .FALSE.) 

CALL  VCURVEdZ,  .FALSE.) 

CALL  WCURVEdZ,  .FALSE.) 
NCART<1> 
XWARNINGSI I  I  I  I  I 


A)  WHEN  USING  BFCS  ST0VAR(H3),  ST0VAR(C4),  ST0VAR(21)  ARE 
AVAILABLE  ONLY  FOR  STORING  NON-ORTHOGONAL  VELOCITY 
COMPONENTS. 

MULTI-RUNS  ARE  NOT  ALLOWED  WITH  BFC  OPTION. 
MOVING  GRID, TWO-PHASE  AND  PARABOLIC  OPTIONS  ARE  NOT 
AVAILABLE  WITH  BFC  OPTION. 

KE-EP  TURBULENCE  MODEL  SHOULD  BE  USED  WITH  BFCS  ONLY 
WHEN  THE  MAIN  FLOW  IS  IN  THE  IZ  DIRECTION. 
BUILT-IN  GRAVITY  TERMS  DO  NOT  TAKE  ACCOUNT  OF  BFCS. 


B) 
C) 

D) 

E) 

XNOTES 


A)  THE  STANDARD  VELOCITY-FIELD  PRINTOUT  FOR  THE 
VELOCITY  RESOLUTES  IS  ACTIVATED  IN  THE  USUAL 
WAY.  AN  ADDITIONAL  OPTION  EXISTS  FOR  PRINTING  THE 
CARTESIAN  VELOCITY-COMPONENTS  WHICH  MAY  BE 
ACTIVATED  BY  SETTING  THE  FOLLOWING  LOGICALS: 

ST0VAR(U2)=.T.  FOR  U-COMPONENT  (CARTESIAN) 
ST0VAR(V2)=.T.  FOR  V-COMPONENT  (CARTESIAN) 
ST0VAR(W2)=.T.  FOR  W-COMPONENT  (CARTESIAN) 

SIMILARLY  PRINTOUT  OF  NON-ORTHOGONAL  VELOCITY 

COMPONENTS  MAY  BE  ACTIVATED  AS  FOLLOWS: 

ST0VAR(C4)=.T.  FOR  U-COMPONENT  (NON-ORTHOG) 
ST0VAR(H3)=.T.  FOR  V-COMPONENT  (NON-ORTHOG) 
ST0VAR(21)=.T.  FOR  W-COMPONENT  (NON-ORTHOG) 

B)  BFC  (TO  ACTIVATE  THE  BFC  OPTION),  IGEN  (THE  CODE  FOR  METHOD 
OF  GRID  SPECIFICATION),  ND  (NUMBER  OF  SUB-DOMAINS)  AND 
NBFC  (THE  Fl  ARRAY  DIMENSION),  MUST  BE  SET  BEFORE 
"STANDARD  BFC  SECTION  2".  ====== 


VAN01450 
VAN01460 
VAN01470 
VAN01480 
VAN01490 
VAN01500 
VAN01510 
VAN01520 
VAN01530 
VAN01540 
VAN01550 
VAN01560 
VAN01570 
VAN01580 
VAN01590 
VAN0160D 
VAN01610 
VAN01620 
VAN0I630 
VAN01640 
VAN01650 
VAN01660 
VAN01670 
VAN01680 
VAN01690 
VAN01700 
VAN01710 
VAN01720 
VAN01730 
VAN01740 
VAN01750 
VAN01760 
VAN01770 
VAN01780 
VAN01790 
VAN0I800 
VAN018I0 
VAN01820 
VAN01830 
VAN01840 
VAN01850 
VAN01860 
VAN01870 
VAN01880 
VAN01890 
VAN01900 
VAN01910 
VAN01920 
VAN01930 
VAN01940 
VAN01950 
VAN01960 
VAN01970 
VAN01980 
VAN01990 
VAN02000 
VAN02010 
VAN02020 
VAN02030 
VAN02040 
VAN02050 
VAN02060 
VAN02070 
VAN02080 
VAN02090 
VAN02100 
VAN02110 
VAN02120 
VAN02130 
VAN02140 
VAN02150 
VAN02I60 
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ALL  OTHER  BFC  DATA  MUST  BE  SET  AFTER  "STANDARD  BFC 
SECTION  2.  ===== 

C)  NXP1,  NYP1,  NZP1  STORE  NX+1,  NY+1,  NZ+1;  THESE  ARE 
AVAILABLE  TO  USER  AFTER  STANDARD  BFC  SECTION  2. 

D)  FOR  IGEN=1  USE  BXFRAC, BYFRAC  &  BZFRAC  IN  PLACE  OF 
XFRACYFRAC  &  ZFRAC. 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  BFC  SECTION  1  STARTS: 
C     DEFAULT  SETTINGSt 

NCART=10 

BFC=.TRUE. 

IGEN=1 

ND=1 

NBFC=5000 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  BFC  SECTION  1  ENDS.. 
C    XUSER  SETS  BFC,  IGEN,  ND  AND  NBFC  HERE: 
C 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  BFC  SECTION  2  STARTS: 

CAL L  SB4I ( NXP1 , NX+1 , NYP1 , NY+1 , NZP1 , NZ+1 , 1 , 0 ) 

IF(BFC)  CALL  BFCDFTCNBFC, XE, XW, YN, YS, ZH,ZL, ND, NXP1, NYP1 , 
&  NZP1,NZ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  BFC  SECTION  2  ENDS. 
C     *USER  SETS  ALL  OTHER  BFC  VARIABLES  HERE: 
C     XUSING  NONIFORM  GRID  1-8 

GTH=65.E-3 

GTL=150.E-3 

GBETA=4. 

GBETA=GBETA*3. 1415927/180 

GTAB=TAN(GBETA) 

DELMAX=2.E-3 

GNBL=5. 

GPWR=2. 

DO  64  IY=1,5 

64  BYFRAC(IY)=(FLOAT(IY)/GNBL)xxGPWR*DELMAX/GTH 
BYFRAC(6)=BYFRAC(5)+3.E-3/GTH 
DEL=(1.-BYFRAC(6))/(FL0AT(NY)-GNBL-1) 

DO  65  IY=7,NY 

65  BYFRAC(IY)=BYFRAC(IY-1)+DEL 

BZFRAC(l)=10.E-3 
DO  66  IZ=2,5 

66  BZFRAC(IZ)=10.E-3+BZFRAC(IZ-l) 
BZFRAC(6)=BZFRAC(5)+5.E-3 

DO  67  IZ=7,9 

67  BZFRAC(IZ)=BZFRAC(IZ-l)+2.E-3 
DO  68  IZ=10,10 

68  BZFRAC(IZ)=BZFRAC(IZ-l)+l.E-3 
DO  77  IZ=11,14 

77  BZFRAC(IZ)=BZFRAC(IZ-l)+.5E-3 
DO  78  IZ=15,15 

78  BZFRAC(IZ)=BZFRAC(IZ-l)+l.E-3 
BZFRAC(16)=BZFRAC(15)+2.E-3 
BZFRAC(17)=BZFRAC(16)+3.E-3 
BZFRAC(18)=BZFRAC(17)+5.E-3 
DO  69  IZ=19,22 

69  BZFRAC(IZ)=BZFRAC(IZ-l)+10.E-3 
BZFRAC(23)=BZFRAC(22)+3.E-3 
BZFRAC( 24 ) =BZFRAC( 23 ) +2 . E-3 
BZFRAC(25)=BZFRACC24)+2.E-3 
BZFRAC(26)=BZFRAC(25)+3.E-3 
BZFRAC( 27 ) =BZFRAC( 26 )+5 . E-3 

DO  71  IZ=28,NZ 

71  BZFRACC  IZ)  =  BZFRAC(  IZ-D  +  10  .  E-3 
DO  72  IZ=1,NZ 

72  BZFRAC(IZ)=BZFRAC(IZ)/GTL 
CALL  DOMAIN(l,l,NX,l,NY,l,NZ) 
DO  61  IX=1,NXP1 

DO  62  IY=1,NYP1 
ZL(IX,IY,1)=0.0 
62  ZH(IX,IY,1)=GTL 
DO  63  IZ=1,NZP1 
YN(IX,IZ,1)=GTH 
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63  YS(IX,IZ,1)=0.0 
C  YSCIX,13,1)  SHOULD  COME  AFTER 
DO  662  IZ=16,25 

YSCIX?izfl)=CBZFRAC(IZ-l)-BZFRAC(3))*GTABxGTL 

DO  663  IZ=13,15 

GZ12=(BZFRAC(IZ-1)-BZFRAC(11))XGTL 
663  YS(IX,IZ,1)=SQRT(YSCIX,16,1)*GZ12*2.-GZ12X*2) 

DO  664  IZ=26,NZ 
YS(IX,IZ,1)=YS(IX,25,1) 

CONTINUE 

STORSA(IFIX(LOW))=.TRUE. 

STORSA(IFIX(HIGH))=.TRUE. 

STORSA(IFIX(SOUTH))=.TRUE. 

STORWD( I FIX( SOUTH) )=. TRUE. 

STORP=.TRUE. 

PRTBFC=.TRUE. 

DARCY=1.E10 


CCC 
662 


664 
61 


CDAR 

C 

C—   GROUP  7 


BLOCKAGE:  BLOCK< . F. >, IPLANE, IPWRIT 


XSET  CONSTANT  POROSITIES  OVER  SUB-DOMAINS  USING: 
CALL  CONPORdR, TYPE, VALUE, IXF, IXL , IYF, IYL , IZF, IZL ) ,  WHERE: 
IR  =  RUN  SECTION  NUMBER,  E.G.  1  FOR  RUN1  SECTION;  'TYPE'=  EAST, 
WEST,  NORTH,  SOUTH,  HIGH,  LOW  &  CELL.  'VALUE' =WANTED  POROSITY 
OVER  REGION  IXF, . . .IZL. 

XDIMENSION  ARRAYS  PE(NX, NY, NZ) ,  PNC  NX, NY, NZ) ,  PH(NX, NY, NZ) ,  & 
PC(NX,NY,NZ)  ABOVE. 

XFOR  FULLY-BLOCKED  CELLS  (IE.  'VALUE»=  0.0)  USER  NEED  SET  ONLY 
THE  'CELL'  POROSITY  (TO  ZERO),  AS  CELL-FACE  AREAS  ARE  THEN 
AUTOMATICALLY  ZEROED. 

XFOR  SATELLITE  PRINTOUT  OF  ALL  POROSITIES  IN  DOMAIN,  »IPLANE'= 
XPLANE  YPLANE  OR  ZPLANE,  FOR  DESIRED  CROSS-SECTION  DIRECTION. 

XFOR  EACH  'TYPE'  A  MAXIMUM  OF  10  CALLS  TO  CONPOR  IS  ALLOWED, 
BUT  IF  REQUIREMENTS  EXCEED  THIS  PROVISION  SET  BLOCK=.T.  & 
IPWRIT=-1,  AND  SET  POROSITY  ARRAYS  EXPLICITLY  HERE  AS  WANTED. 
IN  THIS  CASE,  THE  USER   MUST    SET   A  L  L   ELEMENTS  OF 
ARRAYS  PE,  PN,  PH,  PC  (MANY  MAY  BE  0.0  OR  1.0).  HE  MAY  USE: 
CALL  CR(PARRAY,VALUE,IXF,IXL,IYF,IYL,IZF,IZL,NX,NY,NZ) 
ANY  NUMBER  OF  TIMES,  TO  SET  'PARRAY'  (=  PE,  ETC.)  TO 
•VALUE'  OVER  RANGE  IXF  TO  IXL,  IYF  TO  IYL,  IZF  TO  IZL. 

XCONPOR   MUST   NOT   BE  USED  IN  CONJUNCTION  WITH  EXPLICIT 
SETTINGS  OF  THE  ARRAYS  (INCLUDING  SETTINGS  VIA  CR). 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C—   GROUP  8. DEPENDENT  VARIABLES  TO  BE  SOLVED  FOR  OR  STORED  » 

C      SOLVAR( 1-25 )<25X.F.>,ST0VAR( 1-25 )<25X.F.>,C0NC1( 1-4 )<4X.T.> 

C      USE  FOLLOWING  NAMED  INTEGERS  FOR  ARRAY  ELEMENTS  1-20: 

C      P1,PP,U1,U2,V1,V2,W1,W2,M1,M2,RS,KE,EP,H1,H2,H3,C1,C2,C3,C4. 

S0LVAR(P1)=.TRUE. 

SOLVAR(PP)=.TRUE. 

S0LVAR(V1)=.TRUE. 

S0LVAR(W1)=.TRUE. 

S0LVAR(H1)=.TRUE. 

SOLVAR(KE)=.TRUE. 

SOLVAR(EP)=.TRUE. 

ST0VAR(V2)=.TRUE. 

ST0VAR(W2)=.TRUE. 

ST0VAR(C1)=.TRUE. 

ST0VAR(C2)=.TRUE. 

ST0VAR(C3)=.TRUE. 


C GROUP   9.    VARIABLE   LABELS    : 

C  TITLE(1-25X2HP1,2HPP,2HU1,2HU2,2HV1,2HV2,2HW1,2HW2,2HR1, 

C  2HR2,2HRS,2HKE,2HEP,2HH1,2HH2,2HH3,2HC1,2HC2, 

C      •      2HC3,2HC4,2HRX,2HRY,2HRZ,    2x«Hxxxx> 

TITLE(C1)=TITC1 

TITLE(C2)=TITC2 

TITLE(C3)=TITC3 

TITLE(PP)=TITPP 


C GROUP  10  PROPERTIES: 

C      IRHOK1>,IRH02<1>,RHOK1.0>,RH02<1.0>, 
C     ARHOK1 .  0>,  BRHOK1 .  0>,CRH01<1 .  0> 
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C  IEMUK1>»EMU1<1.0>,EMULAM<1.E-10> 

C  IHSAT,H1SAT,H2SAT,PSATEX<1.0> 

C  SIGMA(1-25X1.0,2.0,1.,1.E10,1.,1.E10,1.,1.E10, 

C  4X1. 0,1. 314, 1.0,1.E10, 10X1. 0> 

IRH01=-1 

PT0T=55.E5 

T0T=555.55 

RAIR=287. 

GAMA=1.35 

CP=RAIR/(1-1/GAMA) 

TW=323. 

HWALL=TWXCP 

HTOT=CPXTOT 

RHTOT=PTOT/TOT/RAIR 

L0GIC(87)=.TRUE. 

ARH01=RHT0T/PT0TXX(1/GAMA) 

BRH01=1./GAMA 
C   TURBULENT  OR  LAMINAR 

IEMU1=2 
C      IEMU1=-1 

JEMU1=IEMU1 

EMUl=l.E-5 

EMULAM=EMU1 

GEMU1=EMU1 

GPR=.7 

SIGMA(24)=GPR 

SIGMA(14)=.9 


C—      GROUP    11    INTER-PHASE   TRANSFER    PROCESSES    . 

C              ICFIP,CFIPS,IMD0T,CMD0T,CA1K1.E6>,CA2K1.E6> 
c 

C—   GROUP  12  SPECIAL  SOURCES  i 

C      ISPCSOC 1-25), AGRAVX,AGRAVY,AGRAVZ,ABUOY, HREF 

c 

C—      GROUP    13    INITIAL    FIELDS    : 
C  FIINIT(1-25X25X1.E-10> 

C  MACH    NO.    OF    FREE   STREAM 

GMACH=3.2 

A=1+(GAMA-1)/2XGMACHXX2 

TE=TOT/A 

RHE=RHT0T/AXX(1/(GAMA-1)) 

PSTAT  =  PT0T/AXX(GAMA/(GAMA-D) 

RH01=ARH01XPSTATXXBRH01 

SONIC=SQRT(GAMAXRAIRXTE) 

WIN=SONICXGMACH 

RKEIN=0.01XWINXX2 

EPIN=0.16xRKEINxxl.5/GTH/2. 

FIINIT(W1)=WIN 

FIINIT(P1)=PSTAT 

FIINIT(H1)=HT0T 

FIINIT(KE)=RKEIN 

FIINIT(EP)=EPIN 


C 

C— 
C 
C 
C 

14 

C 

C— 
C— 

c 
c 
c 


GROUP  14  BOUNDARY/INTERNAL  CONDITIONS  i 
ILOOP1,ILOOPN,XCYCLE<.F.>,PBAR,REGION(1-10X10X.T.> 
XN.B.  ALL  10  REGIONS  ARE  DEFAULTED  .TRUE..   THE  USER  SHOULD 
SET  REGION(I)=. FALSE.  FOR  UNUSED  REGIONS  'I'. 
DO  14  1=1,10 
REGION(I)=. FALSE. 


GROUP  15  TO  24;  REGIONS  1  TO  10 

ONLY  THOSE  REGIONS  ARE  ACTIVE  WHICH  ARE  SPECIFIED  BY  THE 

USER,  PREFERABLY  BY  WAY  OF.- 

CALL  PLACE(IREGN,TYPE,IXF,IXL,IYF,IYL,IZF,IZL)  & 

CALL  COVAL(IREGN,VARBLE,COEFF, VALUE) 

CALL  PLACE(l,LOW,l,NX,l,NY,l,l) 
C  CALL  C0VAL(1,M1,FIXFLU,WINXRHE) 
CDAR   CALL  COVALC1 ,M1 , 1 . E-20, 1 . E+20XWINXRHE) 

GCM=2XGAMA/WIN/(GAMA-1) 

GVM=PTOTXRHE/RHTOT 

CALL  C0VAL(1,M1,GCM,GVM) 

CALL  COVAL(l,Wl,ONLYMS,WIN) 
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C- 
C- 
C 

c 

c- 
c- 
c 

c- 

c- 

c 

c 

c 

c 

c 

c 

c 

c 


c- 

c- 

c 

c 

c 

c- 
c- 

c- 
c- 
c 


CALL  COVAL(l,Hl,ONLYMS,HTOT) 

CALL  COVAL(l,KE,ONLYMS,RKEIN) 

CALL  COVALCl,EP,ONLYMS,EPIN) 

CALL  PLACE(2,HIGH,1,NX,1,NY,NZ,NZ) 

CALL  COVAL(2,M1,FIXVAL,PSTATXO.) 

CALL  COVAL(2,M1,1000XWINXRHE/PSTAT,PSTAT) 

CALL  C0VAL(2,H1,0NLYMS,HT0T) 

E  VANE   IZ(11,NZ) 

5xBYFRAC(l)xGTH) 

1XGTH 

(0.5XDY1) 

(0.5XDY1XSIGMAC24)) 

3, SOUTH, 1, NX, 1,1, 12, NZ) 
CALL  COVAL(3,W1,GOEFF,0.) 
CALL  C0VAL(3,H1,G0EFH,HWALL) 
CALL  COVAL(3,W1,WALL,0.) 
CALL  C0VAL(3,H1,WALL,HWALLX- 
CALL  COVAL(3,KE,WALL,0.) 
CALL  COVAL(3,EP,WALL,0.) 


WALL  ALONG  TH 
GCM=EMUl/(. 
DY1=BYFRAC( 
GOEFF=EMUl/ 
GOEFH=EMUl/ 
CALL  PLACE( 


GROUP  25  GROUND  STATION  t 
GROSTA< . F . >, NAMLST< . F . > 
xNAMLST  ACTIVATES  NAMELIST  IN  GROUND, 
GROSTA=.TRUE. 


GROUP  26  SOLUTION  TYPE  AND  RELATED  PARAMETERS  i 
WHOLEP< . F. >, SUBPST< . F. >, DONACC< . F. > 
WHOLEP=.TRUE. 


GROUP  27  SWEEP  AND  ITERATION  NUMBERS  t 

FSWEEP<1>,LSWEEP<1>,LITHYD<1>,LITC<1>,LITKE<1>,LITH<1>, 

LITER(1-25X9X1,-1,15X1> 

IVELF<1>,NVEL<1>,IVELL<10000>, 

IKEF<1>,NKE<1>,IKEL<10000>, 

IENTF<1>,NENT<1>,IENTL<10000>, 

ICNCF<1>,NCNC<1>,ICNCL<10000>, 

IRH01F<1>,NRHOK1>,IRH01L<10000>, 

IRH02F<1>,NRH02<1>,IRH02L<10000> 

LSWEEP=1201 

GSWP=LSWEEP 

FSWEEP=801 

LITER(PP)=20 

LITER(V1)=5 

LITER(W1)=5 

LITHYD=2 


GROUP  28  TERMINATION  CRITERIA  : 
ENDITC  1-25X9X1.  E-l  0,0.  5, 15X1.  E-10> 
ENDIT(l)=l.E-5 


GROUP  29  RELAXATION   : 

RLXP<1.>,RLXPXY<1.>,RLXPZ<1.>,RLXRH0<1.>,RLXMDT<1.>, 

DTFALSC  3-25X23X1.  El 0> 

DTFALS(Wl)=l.E-5 

DTFALS(Vl)=l.E-5 

DTFALS(KE)=l.E-5 

DTFALS(EP)=l.E-6 

RLXP=.3 


GROUP  30  LIMITS  i 

VELMAX<1.E10>,VELMIN<-1.E10>,RHOMAX<1.E10>,RHOMIN<1.E-10>, 
TKEMAX<1 . E10>,TKEMIN<1 . E-10>, EMUMAX<1 . E10>, EMUMIN<1 . E-10>, 
EPSMAX<1 . E10>, EPSMIN<1 . E-10>, AMDTMX<1 . E10>, AMDTMN<-1 . E10> 
EPSMAX=1.E13 


GROUP  31  SLOWING  DEVICES  «  SL0RH0<1 . >, SL0EMU<1 . > 
SL0RH0=.2 


GROUP  32  PRINT-OUT  OF  VARIABLES  « 
PRINTC1-25X.T.,  .  F. ,  23X  .T  .  >,  SUBWGR<  .  F.  > 
PRINT(C1)=.TRUE. 
PRINT(C2)=.TRUE. 
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PRINT(C3)=.TRUE. 
PRINT(PP)=.TRUE, 


C- 

C- 

C 

C 

C 

C- 

C- 

C 

C 

C- 

C 

C 

C- 

c 

c- 

c- 

c 

c 

c 

c 

c- 

c- 

c 

c 

c- 

c- 

c 

c 


c- 

c- 

c 

c 

c 

c 

c 

c 

c- 

c- 


GROUP  33  MONITOR  PRINT-OUT  i 

IXMON<l>,IYMON<l>,IZMON<l>,NPRMON<l>,NPRMNT<l> 

NPRMON=10 

IYM0N=2 

IZM0N=12 


GROUP  34  FIELD  PRINT-OUT  CONTROL  i 

NPRINT<100>,NTPRIN<100>,NXPRIN<1>,NYPRIN<1>,NZPRIN<1>, 

IZPRF<1>,ISTPRF<1>,IZPRL<10000>,ISTPRL<10000>   • 

NUMCLS<10>,KOUTPT 

NPRINT=LSWEEP 


GROUP  35  TABLE  CONTROL  i 

TABLES<.F.>,NTABLE,NTABVR,LINTAB,NPRTAB,NMON, 

ITAB(l-8),MTABVR(l-8) 


GROUP  36-38  ARE  NOT  DOCUMENTED  IN  THE  INSTRUCTION 
MANUAL  AND  ARE  INTENDED  FOR  MAINTENANCE  PURPOSES  ONLY 
GROUP  36  DEBUG  PRINT-OUT  SLAB  AND  TIME-STEP  t 
IZPRK1>,IZPR2<1>,ISTPRK1>,ISTPR2<1> 


GROUP  37  DEBUG  SWEEP  AND  SUBROUTINES  i 
KEMU,KMAIN,KINDEX,KGEOM,KINPUT,KSODAT,KCOMPF,KSORCE, 
KS0LV1,KS0LV2,KS0LV3,KC0MPP,KADJST,KFLUX,KSHIFT,KDIF, 
KCOMPU, KCOMPV, KCOMPW, KCOMPR, KWALL , KDBRH0<-1>, KDBEXP, KDBMDT 
KDBGEN 


GROUP  38  MONITOR, TEST, AND  FLAG  i 
MONITR< . F . > , FL AG< . F . > , TEST< . T . > , KFL AG<1 > 
END  OF  MAINTENANCE-ONLY  SECTION 


GROUP  39  ERROR  AND  RESIDUAL  PRINT-OUT  : 

IERRP<1000>,RESREF(1,3-24K25*1.>,RESMAP<.F.>, 

RESID(1-25X2X.F.,23X.T.>,K0UTPT 

RESREF(1)=WINXRHE 

RESREF(7)=WINXRESREF(1) 

RESREF(5)=WINXRESREF(1)X0.1 

RESREF(H1)=HT0TXRESREF(1) 

RESREF(KE)=RKEINXRESREF(1) 

RESREF(EP)=EPINXRESREF(1) 

IERRP=LSWEEP/20 

KOUTPT=LSWEEP/20 


GROUP  AO  SPECIAL  DATA  i  LOGICC1 . . 10) , INTGRC1 . . 10) ,RE(21 . . 30), 
NLSP<1>,NISP<1>,NRSP<1>,SPDATA<.F.>,LSPDA(1),ISPDA(1),RSPDA(1) 
USE  FIRST  10  ELEMENTS  OF  ARRAYS  LOGIC  S  INTGR  AND  21ST 
TO  30TH  OF  ARRAY  RE  FOR  TRANSFERRING  SPECIAL  DATA  FROM 
SATELLITE  TO  GROUND,  BUT  IF  REQUIREMENTS  EXCEED  THIS 
PROVISION  SET  SPDATA  =  .T.,  AND  DIMENSION  ARRAYS  LSPDA, 
ISPDA,  RSPDA  ABOVE  AND  IN  GROUND  AS  NEEDED,  AND  SET  HERE 


GROUP  42  RESTARTS  AND  DUMPS 
SAVEM=.TRUE. 
BFPLOT=.TRUE. 
RESTRT=.TRUE. 


SAVEM< . F . > , RESTRT< . F . >, KINPUT 


GROUP  43  GRAFFIC  i 

GRAPHS<  .  F  .  >,  ORTHOG<  .  T  .  >,  ANTSYM,  NPRT<1>,  ITITL<5x<+HXXXX> 

FOR  A  GRAFFIC  RUN,  DIMENSION  PHI1  &  PHI2  AS  FOLLOWS. 

PHIKNXXNYXNZXNM) 

PHI2((NX+2)X(NY+2)X(NZ+2)X(NM+IBLK))  ,  WHERE 

NM=NO.  OF  VARIABLES  STORED  +  DENSITY(-IES) 

IBLK=0  IF  BL0CK=.FALSE.,=4  IF  A  3D  RUN, 

=3  IF  A  2D.YZ  RUN. 


IF(IRUN.EQ.l)  GO  TO  900 
900  CONTINUE 
C™   ALL  RUNS 
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CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  3  ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  <*    STARTS. 


WRITE  GENERAL  DATA  ON  TO  THE  GUSIE1.DTA  TAPE,  ETC... 

IF(SPDATA)  CALL  WRTSPCC LSPDA, NLSP, ISPDA, NISP, RSPDA, NRSP) 
IF(BLOCK)  CALL  WRTPORCPE, PN, PH, PC, NX, NY, NZ, IPLANE) 
IFCBFC)  CALL  WRTBFCC 14, NBFC, XE, XW, YN, YS,ZH,ZL , 
&ND, NX+1 , NY+1 , NZ+1 , NZ, PRTBFC ) 
OLD  PRACTICES  RETAINED  FOR  REFERENCE: 
IF(SPDATA)  CALL  SPCDAT(IRUN) 
IF(BLOCK)  CALL  PORDAT(IRUN) 
IF(GRAPHS)  CALL  SORT(IRUN) 
IF(RESTRT)  GO  TO  902 
DO  901  INDVAR=1,25 

IF(IFIX(FIINIT(INDVAR)+0.1).NE. 10101)  GO  TO  901 
CALL  FLDDAT(IRUN) 
GO  TO  902 
CONTINUE 
CALL  DATAIO(WRT,10) 

IF(MONITR)  CALL  DATAIOCWRT, -6 ) 
CONTINUE 
STOP 
END 
CXXX    IGEN=1  SO  BFCXYZ  NOT  REQUIRED. 
CXXX    COMMENT  OUT  BOTH  VERSIONS. 


901 
902 

999 


SUBROUTINE 

RETURN 

END 


BFCXYZ  (NXP1,NYP1,NZP1) 
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CSDIRECTIVEXXSATLIT    AMI  LEITNER 

C     LAMINAR  SOLUTION  FOR  NWC5    NY=18  NZ=29  YN=GTH 

C     LECSAT  CONVERTED  TO  DIAMSAT 

C  xFILE  NAME:  MODBFCST.FTN 

C  XABSTRACT.  SATELLITE  MODEL  MAIN  PROGRAM.  THIS  VERSION  IS 

C      FOR  USE  WITH  THE  BODY-FITTED  COORDINATE  SCHEME  (SUMMER  1984 

C      VERSION)  PROVIDED  AS  AN  ATTACHMENT  TO  SPRING  1983  PHOENICS. 

C  XDOCUMENTATION:  PHOENICS  INSTRUCTION  MANUAL  (SPRING  1983) 

C     WITH  BODY-FITTED  COORDINATES  INSTRUCTION  SUPPLEMENT 

C      (SUMMER  1984). 

C  ^AUXILIARY  SUBROUTINES  (TAPES,  ETC.)  ARE  IN  SATELLITE  LIBRARY 

C     SERVICEU,  WHICH  MUST  BE  INCLUDED  IN  LINK  EDIT  TO  RUN. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  1  STARTS: 


CHAPTER  1 
C— 


COMMON  BLOCKS  AND  USER'S  DATA. 


1  ENDS. 


INCLUDE  (CMNGUS) 

INCLUDE  (CMNGRF) 
INCLUDE  (GUSSEQ) 
C0MM0N/CPI/IPWRIT,IDUM(243) 
DIMENSION  GDTAPE(3),DFAULT(4) 

DIMENSION  ARRAY1(309),ARRAY2(194),ARRAY3(421) 
LOGICAL  ARRAY1,LSPDA,WRT,RD,NAMLST 
INTEGER  ARRAY2,XPLANE,YPLANE,ZPLANE 

INTEGER  P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS,EP,H1,H2,H3,C1,C2, 
8C3,C4 
REAL  NORTH, LOW 
LOGICAL  BFC 

EQUIVALENCE  ( ARRAY1 ( 1 ) , CARTES ) , ( ARRAY2( 1 ) , NX) 
EQUI VAL ENCE  ( ARRAY3( 1 ) , SPARE1 ( 1 ) ) , (Ml , Rl ) , (M2 , R2) 
EQUIVALENCE  ( LSTRUN, INTGR( 12) ) , (NAMLST, L0GIC(88 ) ) 
EQUIVALENCE  ( LOGIC(20 ) , BFC) 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION 
C$DIRECTIVE*XCMNBF1$$ 

C   THIS  FILE  CONTAINS  SATELLITE  COMMON  BLOCKS  FOR  BFC'S 
C   Fl  MUST  BE  DIMENSIONED  TO  GREATER  THAN  OR  EQUAL  TO 
C  (NX+NY+17XNZ+24XNXXNY+6X(NX+1)X(NY+1)+6*ND) .  THE  VALUE 
C   OF  THE  DIMENSION  MUST  BE  SET  AS  NBFC  IN  GROUP  6  OF  SATLIT. 
COMMON/F0B/FK5000) 
COMMON/CIB/ND/CIC/KOORD 

COMMON/CID/KDBGG,KDBGMF,KDBGCD,KDBIND,KDBMFX,KDBCDT,KDBPCS, 
8  KDBGUV,KDBGPV 

COMMON/CIE/KDBGS,KDBINS 
COMMON/CIF/IGEN/CIG/NCART 
C   THE  FOLLOWING  ARRAYS  MUST  BE  EXACTLY  DIMENSIONED  FOR  NXP1, 
C   NYP1  AND  NZP1,  BUT  MAY  BE  OVER  DIMENSIONED  FOR  ND. 
C   THE  BFRAC  ARRAYS  MUST  BE  DIMENSIONED  TO  ALLOW  FOR  SETTINGS 
C   IN  SATLIT,  THEY  MAY  BE  OVER  DIMENSIONED. 
COMMON/CRA/XW(19,30,1)/CRB/XE(19,30,1) 
&      /CRC/YS(2,30,1)/CRD/YN(2,30,1) 
&      /CRE/ZL(2,19,1)/CRF/ZH(2,19,1) 

&       /CRG/RC0N/CRH/DARCY/CRI/BXFRAC(99)/CRJ/BYFRAC(99) 
&       /CRK/ BZFRAC( 99 ) 

C0MM0N/CLA/ST0RSA(6),ST0RWD(6),ST0RP,ST0RPE,ST0RPN, 
8  ST0RPH,ST0R1,ST0R2,ST0R3,ST0UNV,PRTBFC,ST0CRN 

COMMON/CLC/BFPLOT 

LOGICAL  ST0RP,ST0RPE,ST0RPN,ST0RPH,ST0R1,ST0R2,ST0R3, 
8         STORSA,STORWD,STOUNV,PRTBFC,BFPLOT,STOCRN 
C      END 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  1  STARTS: 
C      GRAFFIC  ARRAYS  DIMENSIONED  AS  NEEDED... 

COMMON/GRAF1/PHIK1)  /GRAF2/PHI2( 1 ) 
C      POROSITY  8  SPECIAL  DATA  ARRAYS  DIMENSIONED  AS  NEEDED... 

DIMENSION  PE(1,1,1),PN(1,1,1),PH(1,1,1),PC( 1,1,1) 

DIMENSION  LSPDA(1),ISPDA(1),RSPDA(1) 
C      USER  PLACES  HIS  VARIABLES,  ARRAYS,  EQUIVALENCES  ETC.  HERE. 

EQUIVALENCE(RAIR,RE(21)),(GAMA,RE(22)),(GSWP,RE(23)) 
1,(GPR,RE(24)),(TW,RE(25)),(GEMU1,RE(26)),(JEMU1,INTGR(1)) 
C      USER  PLACES  HIS  DATA  STATEMENTS  HERE. 

DATA  NLSP,NISP,NRSP/1,1,1/ 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  1  ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  2  STARTS: 
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CHAPTER  2    SET  CONSTANTS,  AND  ARRANGE  FILE  MANIPULATIONS. 

C      PLEASE~DO~NOT  ALTER,  OR  RE-SET,  ANY  OF  THE  REMAINING 
C      STATEMENTS  OF  THIS  CHAPTER. 

DATA  CELL, EAST, WEST, NORTH, SOUTH, HIGH, LOW, VOLUME/ 
&  0.,1.,2.,3.,4.,5.,6. ,7.  / 

DATA  P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS,KE,EP,H1,H2,H3,C1,C2, 
&C3,C4/1, 2, 3,4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20/ 

DATA  FIXFLU,FIXVAL,ONLYMS,WALL/l.E-10,l.E10,0.0,-10.0/ 

DATA  IPLANE,XPLANE,YPLANE,ZPLANE/0,1,2,3/ 

DATA  WRT,RD,DFAULT/.TRUE., . FALSE. , 4HDEFA, 4HULT . ,4HDTA/, 1HG/ 

DATA  GDTAPE/4HGUSI,4HE1.D,2HTA/ 

DATA  NLDATA,NIDATA,NRDATA/309,194,421/ 

DATA  NLCREG,NTCVRG/60,350/ 

DATA  TITPP,TITC1,TITC2,TITC3/3HRH0,4HMACH,4HTEMP,4HCFST/ 

CALL  TAPES(10,GDTAPE,3,1,4*NRDATA) 

c READ  DEFAULT  FILE  IF  BLOCKDATA  ABSENT 

IF(INTGR1(29).NE.10)  GO  TO  2 

CALL  WRIT40C40HDATA  ESTABLISHED  IN  BLOCK  DATA.  ) 

GO  TO  3 
2  CALL  DEFLT 

CALL  TAPES(1,DFAULT,4,2,4XNRDATA) 

CALL  DATAIO(RD,l) 

CALL  WRIT40C40HDATA  TAKEN  FROM  DEFAULT. DTA  ON  GROUP  A/C) 

CALL  WRIT40C40HFILE  MODSTL.FTN  IS  THE  SATLIT  USED.      ) 

L0GIC(89)=.TRUE. 


CD 
CD 


CHAPTER  3    DEFINE  DATA  FOR  NRUN  RUNS. 

c 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  2  ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  2  STARTS: 
C—  GROUP   41MULTI-RUNS  :  RUNC1-30K  .  T  .  ,  29*  .  F.  > 
C 

RUN(1)=. FALSE. 
C   NOTE:  ALL  RUNS  ARE  DEACTIVATED  AT  THIS  POINT  -  USER  SHOULD 
C   ====   SWITCH  ON  ONE  ONLY  OF  RUNS  1-4  IN  NEXT  STATEMENT. 

RUN(1)=.TRUE. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  2  ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  3  STARTS: 
DO  10  IRUN=1,30 

IFC.NOT.RUN(IRUN))  GO  TO  10 

NRUN=NRUN+1 

LSTRUN=IRUN 
10  CONTINUE 

DO  999  IRUN=1,LSTRUN 

IFC.NOT.RUN(IRUN))  GO  TO  999 

INTGR(ll)  =  IRUN 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  3  ENDS. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  3  STARTS: 

C ALL  INTEGER  VARIABLES  ARE  DEFAULTED  TO  0,  AND  REAL  VARIABLES 

C      TO  0.0,  UNLESS  OTHERWISE  INDICATED. 

C      E.G.  BY  VARIABLE<10>,  OR  <10.0>  AS  APPROPRIATE. 

C      THE  DEFAULT  SETTINGS  OF  ALL  LOGICAL  VARIABLES  ARE  ALWAYS 

C      INDICATED,  E.G.  VARIABLE< .T . >,  OR   VARIABLE< . F . > . 

C 

C—   RUN1 

c 

C—   GROUP  1.  FLOW  TYPE  : 

C      PARAB< . F . > , CARTES< . T . > , ONEPHS< . T . > 

c 

C—   GROUP  2.  TRANSIENCE  : 


C 

C 

c 
c 

c 

c— 

c 

c 

c 

c 


STEADY<.T.>,ATIME,LSTEP<1>,FSTEP<1> 

TLAST<1.E10>,TFRAC(1-30K30X1.> 

SERVICE  SUBROUTINE  FOR  'NT1  POWER-LAW  TIME  STEPS: 

CALL  GRDPWRCO, NT, TLAST, POWER) 


GROUP  3.  X-DIRECTION  : 
NX<l>,XULAST<1.0>,XFRAC(l-30) 
SERVICE  SUBROUTINE  FOR  POWER-LAW  GRID: 
CALL  GRDPWRd, NX, XULAST, POWER) 
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C- 
C 
C 
C 

C- 

C- 

C 

C 

C 

C- 
C- 

c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


GROUP  4.  Y-DIRECTION  : 

NY<1>,YVLAST<1.0>,YFRAC(1-30),RINNER,SNALFA 
SERVICE  SUBROUTINE  FOR  POWER-LAW  GRID: 
CALL  GRDPWRC2, NY, YVLAST, POWER) 
NY=18 


GROUP  5.  Z-DIRECTION  i 

NZ<l>,ZWLAST<1.0>,ZFRAC(l-30) 

SERVICE   SUBROUTINE   FOR   POWER-LAW   GRID: 

CALL    GRDPWR(3,NZ,ZWLAST, POWER) 

NZ=29 


0>, 


GROUP   6.    MOVING   GRID   OR   DISTORTED    (BODY-FITTED)    GRID    « 

—   MOVING   GRID    i 
MGRID,IZW1,IZW2,AZW2,BZW2,CZW2,PINT,ZW2M1T 

BODY-FITTED   GRID    — 

BFC<.T.>,IGEN<1>,ND<1>,NBFC<500  0>,KOORD,RCON 

BXFRAC(1-NXK1.0,NXM1X0.0> 

BYFRAC(1-NYX1.0,NYM1X0.0> 

BZFRACC 1-NZX1 .  0,  NZM1*0  .  0> 

SERVICE    SUBROUTINE    FOR    SUB-DOMAIN    SPECIFICATION    (FOR    IGEN=1 

ONLY): 

CALL    DOMAIN(ID,IXF,IXL,IYF,IYL,IZF,IZL) 
XE(1-NYP1,1-NZP1,1-ND)<(NYP1XNZP1XND)*1.0>, 
XW(1-NYP1,1-NZP1,1-ND), 

YN(1-NXP1,1-NZP1,1-NDX(NXP1XNZP1XND)X1 
YS(1-NXP1,1-NZP1,1-ND), 

ZH(1-NXP1,1-NYP1,1-NDX(NXP1XNYP1XND)X1.0>, 
ZL(1-NXP1,1-NYP1,1-ND),ST0RSA(1-6X6X.F.>,ST0RWD(1-6X6X, 
STORP<.F.>,STORPE<.F.>,STORPN<.F.>,STORPH<.F.>,STOUNV<.F.>, 
PRTBFC<.F.>,DARCY,BFPLOT<.F.> 

CYCLIC  BOUNDARY  CONDITIONS  ARE  DEFAULTED  INACTIVE  ; 
TO  ACTIVATE  THEM  AT  SELECTED  IZ  SLABS  USE  SERVICE  SUBROUTINE: 

CALL  XCYIZdZ,  .TRUE.  ) 
SERVICE  SUBROUTINE  TO  DEACTIVATE  CURVATURE  TERMS  IN  U,  V 
AND  W  EQUATIONS  ASSOCIATED  WITH  CURVATURE  OF  IX,  IY,  IZ 
GRID  LINES  RESPECTIVELY: 

CALL  UCURVEdZ,  .FALSE.) 

CALL  VCURVEdZ,  .FALSE.) 

CALL  WCURVEdZ,  .FALSE.) 
NCART<1> 
XWARNINGSI I | |  |  | 


F.>, 


A)  WHEN  USING  BFCS  ST0VAR(H3),  ST0VAR(C4),  ST0VAR(21)  ARE 
AVAILABLE  ONLY  FOR  STORING  NON-ORTHOGONAL  VELOCITY 
COMPONENTS. 

MULTI-RUNS  ARE  NOT  ALLOWED  WITH  BFC  OPTION. 
MOVING  GRID, TWO-PHASE  AND  PARABOLIC  OPTIONS  ARE  NOT 
AVAILABLE  WITH  BFC  OPTION. 

KE-EP  TURBULENCE  MODEL  SHOULD  BE  USED  WITH  BFCS  ONLY 
WHEN  THE  MAIN  FLOW  IS  IN  THE  IZ  DIRECTION. 
BUILT-IN  GRAVITY  TERMS  DO  NOT  TAKE  ACCOUNT  OF  BFC'S. 


B) 
C) 

D) 

E) 
XNOTES 


A)  THE  STANDARD  VELOCITYcFIELD  PRINTOUT  FOR  THE 
VELOCITY  RESOLUTES  IS  ACTIVATED  IN  THE  USUAL 
WAY.  AN  ADDITIONAL  OPTION  EXISTS  FOR  PRINTING  THE 
CARTESIAN  VELOCITY-COMPONENTS  WHICH  MAY  BE 
ACTIVATED  BY  SETTING  THE  FOLLOWING  LOGICALS: 

ST0VAR(U2)=.T.  FOR  U-COMPONENT  (CARTESIAN) 
ST0VAR(V2)=.T.  FOR  V-COMPONENT  (CARTESIAN) 
ST0VAR(W2)=.T.  FOR  W-COMPONENT  (CARTESIAN) 

SIMILARLY  PRINTOUT  OF  NON-ORTHOGONAL  VELOCITY 

COMPONENTS  MAY  BE  ACTIVATED  AS  FOLLOWS: 

ST0VAR(C4)=.T.  FOR  U-COMPONENT  (NON-ORTHOG) 
ST0VAR(H3)=.T.  FOR  V-COMPONENT  (NON-ORTHOG) 
ST0VAR(21)=.T.  FOR  W-COMPONENT  (NON-ORTHOG) 

B)  BFC  (TO  ACTIVATE  THE  BFC  OPTION),  IGEN  (THE  CODE  FOR  METHOD 
OF  GRID  SPECIFICATION),  ND  (NUMBER  OF  SUB-DOMAINS)  AND 
NBFC  (THE  Fl  ARRAY  DIMENSION),  MUST  BE  SET  BEFORE 
"STANDARD  BFC  SECTION  2".  ====== 
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ALL  OTHER  BFC  DATA  MUST  BE  SET  AFTER  "STANDARD  BFC 
SECTION  2.  ===== 

C)  NXP1,  NYP1,  NZP1  STORE  NX+1,  NY+1,  NZ+1;  THESE  ARE 
AVAILABLE  TO  USER  AFTER  STANDARD  BFC  SECTION  2. 

D)  FOR  IGEN=1  USE  BXFRAC, BYFRAC  &  BZFRAC  IN  PLACE  OF 
XFRACYFRAC  &  ZFRAC. 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  BFC  SECTION  1  STARTS 
C      DEFAULT  SETTINGS. 

NCART=10 

BFC=.TRUE. 

IGEN=1 

ND  =  1 

NBFC=5000 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  BFC  SECTION  1  ENDS.: 
C     XUSER  SETS  BFC,  IGEN,  ND  AND  NBFC  HERE: 
C 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  BFC  SECTION  2  STARTS: 

CALL  SB4I  ( NXP1 , NX+1 , NYP1 , NY+1 , NZP1 , NZ+1 ,1,0) 

IF(BFC)  CALL  BFCDFT(NBFC,XE, XW, YN, YS,ZH,ZL , ND, NXP1 , NYP1 , 
&  NZP1,NZ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  BFC  SECTION  2  ENDS. 
C     XUSER  SETS  ALL  OTHER  BFC  VARIABLES  HERE: 
C     XUSING  NONIFORM  GRID  1-8 

GTH=65.E-3 

GTL=150.E-3 

GBETA=4. 

GBETA=GBETA*3. 1415927/180 

GTAB=TAN(GBETA) 

DELMAX=2.E-3 

GNBL=5. 

GPWR=4. 

DO  64  IY=1,5 

64  BYFRAC(IY)=(FLOAT(IY)/GNBL)*XGPWR*DELMAX/GTH 
BYFRAC(6)=BYFRAC(5)+3.E-3/GTH 
DEL=(1.-BYFRAC(6))/(FL0AT(NY)-GNBL-1) 

DO  65  IY=7,NY 

65  BYFRAC(IY)=BYFRAC(IY-1)+DEL 

BZFRAC(l)=10.E-3 
DO  66  IZ=2,5 

66  BZFRAC(IZ)=10.E-3+BZFRAC(IZ-l) 
BZFRACC6)=BZFRACC5)+5.E-3 

DO  67  IZ=7,9 

67  BZFRAC(IZ)=BZFRAC(IZ-l)+2.E-3 
DO  68  IZ=10,11 

68  BZFRAC(IZ)=BZPRAC(IZ-l)+.5E-3 
BZFRAC(12)=BZFRAC(ll)+l.E-3 
DO  77  IZ=13,14 

77  BZFRAC(IZ)=BZFRAC(IZ-l)+.5E-3 
DO  78  IZ=15,15 

78  BZFRAC(IZ)=BZFRAC(IZ-l)+l.E-3 
BZFRAC(16)=BZFRAC(15)+l.E-3 
BZFRAC(17)=BZFRAC(16)+2.E-3 
BZFRAC(18)=BZFRAC(17)+7.E-3 
DO  69  IZ=19,22 

69  BZFRAC(IZ)=BZFRAC(IZ-l)+10.E-3 
BZFRAC(23)=BZFRAC(22)+3.E-3 
BZFRAC(24)=BZFRAC(23)+2.E-3 
BZFRAC(25)=BZFRAC(24)+2.E-3 
BZFRAC(26)=BZFRAC(25)+3.E-3 
BZFRAC(27)=BZFRAC(26)+5.E-3 

DO  71  IZ=28,NZ 

71  BZFRAC(IZ)=BZFRAC(IZ-l)+10.E-3 
DO  72  IZ=1,NZ 

72  BZFRAC(IZ)=BZFRAC(IZ)/GTL 
CALL  D0MAIN(1,1,NX,1,NY,1,NZ) 
DO  61  IX=1,NXP1 

DO  62  IY=1,NYP1 

ZL(IX,IY,1)=0.0 

62   ZH(IX,IY,1)=GTL 

DO  63  IZ=1,NZP1 
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YN(IX,IZ,1)=GTH 
63   YS(IX,IZ,1)=0.0 
C   YS(IX,13,1)  SHOULD  COME  AFTER 

DO  662  IZ=5,25 
CBL    DO  662  IZ=16,25 

662   YS(IX,IZ,1)=(BZFRAC(IZ-1)-BZFRAC(3))XGTABXGTL 
CBL    DO  663  IZ=13,15 

CBL    GZ12=(BZFRACCIZ-l)-BZFRACCll))*GTL-.5E-3 
CBL663YS(IX,IZ,l)=SQRT(YS(IX,16,l)xGZ12x2.-GZ12xx2) 
DO  664  IZ=26,NZ 
664   YS(IX,IZ,1)=YS(IX,25,1) 
61   CONTINUE 

STORSA(IFIX(LOW))=.TRUE. 

STORSA( IFIX( HIGH) )=. TRUE. 

STORSA(IFIX( SOUTH ))=. TRUE. 

STORWD(IFIX(SOUTH))=.TRUE.   -' 

STORP=.TRUE. 

PRTBFC=.TRUE. 

DARCY=1.E10 


CDAR 

C 

C 

C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 

c 

c 
c 
c 


CT 
CT 


GROUP  7.  BLOCKAGE:  BLOCK< . F. >, IPLANE, IPWRIT 

XSET  CONSTANT  POROSITIES  OVER  SUB-DOMAINS  USINGt 
CALL  CONPORdR,  TYPE,  VALUE,  IXF, IXL , IYF, IYL , IZF, IZL ) ,  WHERE: 
IR=RUN  SECTION  NUMBER,  E.G.  1  FOR  RUN1  SECTION;  «TYPE'=  EAST, 
WEST,  NORTH,  SOUTH,  HIGH,  LOW  &  CELL.  'VALUE' =WANTED  POROSITY 
OVER  REGION  IXF, . . .IZL. 

XDIMENSION  ARRAYS  PECNX, NY, NZ) ,  PN(NX, NY, NZ) ,  PH(NX, NY, NZ) ,  & 
PC(NX,NY,NZ)  ABOVE. 

XFOR  FULLY-BLOCKED  CELLS' (IE.  »VALUE'=  O.O)  USER  NEED  SET  ONLY 
THE  'CELL'  POROSITY  (TO  ZERO),  AS  CELL-FACE  AREAS  ARE  THEN 
AUTOMATICALLY  ZEROED. 

XFOR  SATELLITE  PRINTOUT  OF  ALL  POROSITIES  IN  DOMAIN,  'IPLANE'= 
XPLANE  YPLANE  OR  ZPLANE,  FOR  DESIRED  CROSS-SECTION  DIRECTION. 

XFOR  EACH  'TYPE'  A  MAXIMUM  OF  10  CALLS  TO  CONPOR  IS  ALLOWED, 
BUT  IF  REQUIREMENTS  EXCEED  THIS  PROVISION  SET  BLOCK=.T.  & 
IPWRIT=-1,  AND  SET  POROSITY  ARRAYS  EXPLICITLY  HERE  AS  WANTED. 
IN  THIS  CASE,  THE  USER   MUST    SET   A  L  L   ELEMENTS  OF 
ARRAYS  PE,  PN,  PH,  PC  (MANY  MAY  BE  0.0  OR  1.0).  HE  MAY  USE: 
CALL  CR(PARRAY,VALUE,IXF,IXL,IYF,IYL,IZF,IZL,NX,NY,NZ) 
ANY  NUMBER  OF  TIMES,  TO  SET  'PARRAY'  (=  PE,  ETC.)  TO 
'VALUE'  OVER  RANGE  IXF  TO  IXL,  IYF  TO  IYL,  IZF  TO  IZL. 

XCONPOR   MUST   NOT   BE  USED  IN  CONJUNCTION  WITH  EXPLICIT 
SETTINGS  OF  THE  ARRAYS  (INCLUDING  SETTINGS  VIA  CR) . 


GROUP  8. DEPENDENT  VARIABLES  TO  BE  SOLVED  FOR  OR  STORED  : 

S0LVAR(1-25K25X.F.>,ST0VAR(1-25X25X.F.>,C0NC1(1-<4)<4X.T.> 

USE  FOLLOWING  NAMED  INTEGERS  FOR  ARRAY  ELEMENTS  1-20: 

P1,PP,U1,U2,V1,V2,W1,W2,M1,M2,RS,KE,EP,H1,H2,H3,C1,C2,C3,C4, 

S0LVAR(P1)=.TRUE. 

SOLVAR(PP)=.TRUE. 

S0LVAR(V1)=.TRUE. 

S0LVAR(W1)=.TRUE. 

S0LVAR(H1)=.TRUE. 

SOLVAR(KE)=.TRUE. 

SOLVAR(EP)=.TRUE. 

ST0VAR(V2)=.TRUE. 

ST0VAR(W2)=.TRUE. 

ST0VAR(C1)=.TRUE. 

ST0VAR(C2)=.TRUE. 

ST0VAR(C3)=.TRUE. 


C- 

c- 
c 


GROUP  9.  VARIABLE  LABELS  : 

TITLE(1-25X2HP1,2HPP,2HU1,2HU2,2HV1,2HV2,2HW1,2HW2,2HR1, 

2HR2 , 2HRS, 2HKE, 2HEP, 2HH1 , 2HH2, 2HH3, 2HC1 , 2HC2, 

2HC3, 2HC4, 2HRX, 2HRY, 2HRZ,  2X4Hxxxx> 

TITLE(C1)=TITC1 

TITLE(C2)=TITC2 

TITLE(C3)=TITC3 

TITLE(PP)=TITPP 


GROUP  10  PROPERTIES: 
IRH0K1>,IRH02<1>,RH0K1.0>,RH02<1.0>, 
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r  ARHOK1.0>,BRH01<1.0>,CRH01<1.0> 

r  IEMUK1>,EMUK1.0>,EMULAM<1.E-10> 

C     IHSAT,H1SAT,H2SAT,PSATEX<1.0> 

C     SIGMAC  1-25X1.0,  2.  0,1.,1.  El  0, 1 . ,  1 . E10, 1 . , 1 . E10, 

C     4X1. 0,1. 314, 1.0,1.E10, 10X1. 0> 

IRH01=-1 

PT0T=55.E5 

T0T=555.55 

RAIR=287. 

GAMA=1.35 

CP=RAIR/(1-1/GAMA) 

TW=323. 

HWALL=TWXCP 

HTOT=CPXTOT 

RHTOT=PTOT/TOT/RAIR 

L0GIC(87)=.TRUE. 

ARH01=RHT0T/PT0TXX(1/GAMA) 

BRH01=1./GAMA 
C   TURBULENT  OR  LAMINAR 

IEMU1=-1 
C      IEMU1=1 

JEMU1=IEMU1 

EMUl=l.E-5 

EMULAM=EMU1 

GEMU1=EMU1 

GPR=.7 

SIGMAC24XGPR 

SIGMA(14)=GPR 


C—      GROUP    11    INTER-PHASE   TRANSFER    PROCESSES    : 

C             ICFIP,CFIPS,IMD0T,CMD0T,CA1K1.E6>,CA2K1.E6> 
c 

C—   GROUP  12  SPECIAL  SOURCES  : 

C      ISPCSOC 1-25 ),AGRAVX,AGRAVY,AGRAVZ,ABUOY, HREF 
c 

C—   GROUP  13  INITIAL  FIELDS  : 
C      FIINIT(1-25X25X1.E-10> 
C     MACH  NO.  OF  FREE  STREAM 

GMACH=3.2 

A=1+(GAMA-1)/2XGMACHXX2 

TE=TOT/A 

RHE  =  RHT0T/AXX(1/(GAMA-D) 

PSTAT  =  PTOT/AXX(GAMA/(GAMA-D) 

RH01=ARH01XPSTATXXBRH01 

SONIC=SQRT(GAMAXRAIRXTE) 

WIN=SONICXGMACH 

RKEIN=0.01XWINXX2 

EPIN=0.16xRKEINXxi.5/GTH/2. 

FIINIT(W1)=WIN 

FIINIT(P1)=PSTAT 

FIINIT(H1)=HT0T 

FIINIT(KE)=RKEIN 

FIINIT(EP)=EPIN 


C 

C— 

c 
c 
c 

14 

C 

C— 
C— ^ 
C 

c 
c 

c 

CDAR 


GROUP  14  BOUNDARY/INTERNAL  CONDITIONS  : 
ILOOP1,ILOOPN,XCYCLE<.F.>,PBAR,REGION(1-10X10X.T.> 
XN.B.  ALL  10  REGIONS  ARE  DEFAULTED  .TRUE..   THE  USER  SHOULD 
SET  REGION(I)=. FALSE.  FOR  UNUSED  REGIONS  "I*. 
DO  14  1=1,10 
REGIONCIX.  FALSE. 


GROUP  15  TO  24;  REGIONS  1  TO  10 

ONLY  THOSE  REGIONS  ARE  ACTIVE  WHICH  ARE  SPECIFIED  BY  THE 

USER,  PREFERABLY  BY  WAY  OFt- 

CALL  PLACE(IREGN,TYPE,IXF,IXL,IYF,IYL,IZF,IZL)  & 

CALL  COVAL(IREGN,VARBLE,COEFF, VALUE) 

CALL  PLACE(l,LOW,l,NX,l,NY,l,l) 

CALL  C0VAL(1,M1,FIXFLU,WINXRHE) 

CALL  COVALC 1 ,M1 , 1 . E-20, 1 . E+20XWINXRHE) 

GCM=2XGAMA/WIN/(GAMA-1) 

GVM=PTOTXRHE/RHTOT 

CALL  C0VAL(1,M1,GCM,GVM) 
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CT 
CT 

C- 

c- 

c 

c 


c- 
c- 
c 

c- 

c- 

c 

c 

c 

c 

c 

c 

c 

c 


CR 


CALL  COVAL(l,Wl,ONLYMS,WIN) 

CALL  COVAL(l,Hl,ONLYMS,HTOT) 

CALL  COVAL(l,KE,ONLYMS,RKEIN) 

CALL  COVAL(l,EP,ONLYMS,EPIN) 

CALL  PLACEC2,HIGH,1,NX,1,NY,NZ,NZ) 

CALL  C0VAL(2,M1,FIXVAL,PSTAT*0.) 

CALL  COVAL(2,Ml,1000XWINxRHE/PSTAT,PSTAT) 

CALL  C0VAL(2,H1,0NLYMS,HT0T) 

E  VANE   IZC11,NZ) 

5xBYFRAC(l)xGTH) 

1XGTH 

(0.5XDY1) 

(0.5xDYlxSIGMA(24)) 

3,S0UTH,1,NX,1,1,4,NZ) 
CALL  COVAL(3,W1,GOEFF,0. ) 
CALL  C0VAL(3,H1,G0EFH,HWALL) 
CALL  COVAL(3,W1,WALL,0.) 
CALL  C0VAL(3,H1,WALL,HWALL) 
CALL  COVAL(3,KE,WALL,0.) 
CALL  COVAL(3,EP,WALL,0.) 


WALL  ALONG  TH 
GCM=EMUl/(. 
DY1=BYFRAC( 
G0EFF=EMU1/ 
G0EFH=EMU1/ 
CALL  PLACEC 


GROUP  25  GROUND  STATION  : 
GROSTA< . F . >, NAMLST< . F . > 
XNAMLST  ACTIVATES  NAMELIST  IN  GROUND. 
GROSTA=.TRUE. 


GROUP  26  SOLUTION  TYPE  AND  RELATED  PARAMETERS  i 

WHOLEP<.F.>,SUBPST<.F.>,DONACC<.F.> 

WHOLEP=.TRUE. 


GROUP  27  SWEEP  AND  ITERATION  NUMBERS  : 

FSWEEP<1>,LSWEEP<1>,LITHYD<1>,LITC<1>,LITKE<1>,LITH<1>, 

L  ITER(  1-25X9X1, -1,1 5X1  > 

IVELF<1>,NVEL<1>,IVELL<10000>, 

IKEF<l>,NKE<l>,IKEL<10O0O>, 

IENTF<1>,NENT<1>,IENTL<10000>, 

ICNCF<1>,NCNC<1>,ICNCL<10  000>, 

IRH01F<1>,NRH01<1>,IRH01L<10000>, 

IRH02F<1>,NRH02<1>,IRH02L<10000> 

LSWEEP=400 

GSWP=LSWEEP 

FSWEEP=200 

LITER(PP)=20 

LITER(V1)=5 

LITER(W1)=5 

LITHYD=2 


GROUP  28  TERMINATION  CRITERIA  : 

ENDIT(1-25X9X1.E-10,0.5,15X1.E-10> 

ENDIT(l)=l.E-5 


GROUP    29    RELAXATION       : 

RLXP<1.>,RLXPXY<1.>,RLXPZ<1.>,RLXRH0<1.>,RLXMDT<1.>, 

DTFALS(3-25X23X1.E10> 

DTFALS(Wl)=l.E-5 

DTFALS(Vl)=l.E-5 

RLXP=.2 


GROUP  30  LIMITS  : 

VELMAX<1.E10>,VELMIN<-1.E10>,RHOMAX<1.E10>,RHOMIN<1.E-10>, 
TKEMAX<1.E10>,TKEMIN<1.E-10>,EMUMAX<1.E10>,EMUMIN<1.E-10>, 
EPSMAX<1 . E10>, EPSMIN<1 . E-10>, AMDTMX<1 . E10>, AMDTMN<-1 . E10> 


GROUP  31  SLOWING  DEVICES  i  SL0RH0<1 . >, SL0EMU<1 . > 
SL0RH0=.2 


c 

C—   GROUP  32  PRINT-OUT  OF  VARIABLES 
C 


PRINTU-25X.T., 
PRINT(C1)=.TRUE. 
PRINT(C2)=.TRUE. 
PRINT(C3)=.TRUE. 
PRINT(PP)=.TRUE. 


F.,23X.T.>,SUBWGR<.F.> 
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C- 
C- 
C 
C 

C 

C- 
C- 
C 

c 

c- 

c 

c 

c- 

c 

c- 

c- 

c 

c 

c 

c 

c- 

c- 

c 

c 

c- 

c- 

c 

c 


c- 

c- 

c 

c 

c 

c 

c 

c 


GROUP  33  MONITOR  PRINT-OUT  « 
IXMON<l>,IYMON<l>,IZMON<l>,NPRMON<l>,NPRMNT<l> 

NPRMON=5 

IYM0N=2 

IZM0N=12 


GROUP  34  FIELD  PRINT-OUT  CONTROL  » 

NPRINT<100>,NTPRIN<100>,NXPRIN<1>,NYPRIN<1>,NZPRIN<1>, 

IZPRF<1>,ISTPRF<1>,IZPRL<10000>,ISTPRL<10000> 

NUMCLS<10>,KOUTPT 

NPRINT=LSWEEP 


GROUP  35  TABLE  CONTROL  « 

TABLES<.F.>,NTABLE,NTABVR,LINTAB,NPRTAB,NMON, 

ITAB(l-8),MTABVR(l-8) 


GROUP  36-38  ARE  NOT  DOCUMENTED  IN  THE  INSTRUCTION 
MANUAL  AND  ARE  INTENDED  FOR  MAINTENANCE  PURPOSES  ONLY 
GROUP  36  DEBUG  PRINT-OUT  SLAB  AND  TIME-STEP  « 
IZPRK1>,IZPR2<1>,ISTPRK1>,ISTPR2<1> 


GROUP  37  DEBUG  SWEEP  AND  SUBROUTINES  i 

KEMU, KMAIN, KINDEX, KGEOM, KINPUT, KSODAT, KCOMPF, KSORCE, 

KS0LV1,KS0LV2,KS0LV3,KC0MPP,KADJST,KFLUX,KSHIFT,KDIF, 

KCOMPU, KCOMPV, KCOMPW, KCOMPR, KWALL , KDBRH0<-1>, KDBEXP, KDBMDT 

KDBGEN 


GROUP  38  MONITOR, TEST, AND  FLAG  i 
MONITR< . F . > , FLAG< . F . >, TEST< . T . >,  KFLAG<1> 
END  OF  MAINTENANCE-ONLY  SECTION 


GROUP  39  ERROR  AND  RESIDUAL  PRINT-OUT  i 

IERRP<1000>,RESREF(l,3-24X25xl.>,RESMAP<.F.>, 

RESIDC1-25X2*.F.,23X.T.>,K0UTPT 

RESREF(1)=WINXRHE 

RESREF(7)=WINXRESREFC1) 

RESREF(5)=WIN*RESREF(1)X0.1 

RESREF(H1)=HT0TXRESREF(1) 

RESREF(KE)=RKEINXRESREFU) 

RESREF(EP)=EPINXRESREF(1) 

IERRP=LSWEEP/10 

KOUTPT=LSWEEP/10 


GROUP  40  SPECIAL  DATA  :  LOGIC( 1 . . 10) , INTGRC 1 . . 10 ) , REC21 . . 30) , 
NLSP<1>,NISP<1>,NRSP<1>,SPDATA<.F.>,LSPDA(1),ISPDA(1),RSPDA(1) 
USE  FIRST  10  ELEMENTS  OF  ARRAYS  LOGIC  8  INTGR  AND  21ST 
TO  30TH  OF  ARRAY  RE  FOR  TRANSFERRING  SPECIAL  DATA  FROM 
SATELLITE  TO  GROUND,  BUT  IF  REQUIREMENTS  EXCEED  THIS 
PROVISION  SET  SPDATA  =  .T.,  AND  DIMENSION  ARRAYS  LSPDA, 
ISPDA,  RSPDA  ABOVE  AND  IN  GROUND  AS  NEEDED,  AND  SET  HERE 

SAVEM< . F . > , RESTRT< . F . > , KINPUT 


C—   GROUP  42  RESTARTS  AND  DUMPS 

SAVEM=.TRUE. 

BFPLOT=.TRUE. 
C      RESTRT=.TRUE. 

c :- 

c 

C GROUP  43  GRAFFIC  t 

C  GRAPHS< . F. >, ORTHOG< . T . >, ANTSYM, NPRT<1>, ITITL<5*4HXXXX> 

C—  FOR  A  GRAFFIC  RUN,  DIMENSION  PHI1  8  PHI2  AS  FOLLOWS: 

C  PHIKNXXNYXNZXNM) 

C  PHI2((NX+2)X(NY+2)X(NZ+2)X(NM+IBLK))  ,  WHERE 

C  NM=NO.  OF  VARIABLES  STORED  +  DENSITY(-IES) 

C  IBLK=0  IF  BL0CK=.FALSE.,=4  IF  A  3D  RUN, 

C      =3  IF  A  2D.YZ  RUN. 

c 

IF(IRUN.EQ.l)  GO  TO  900 
900  CONTINUE 
C—   ALL  RUNS 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  3  ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  4  STARTS* 
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C   WRITE  GENERAL  DATA  ON  TO  THE  GUSIE1.DTA  TAPE,  ETC... 

IF(SPDATA)  CALL  WRTSPC( LSPDA, NLSP, ISPDA, NISP, RSPDA, NRSP) 
IF(BLOCK)  CALL  WRTPORCPE, PN, PH, PC, NX, NY, NZ, IPLANE) 
IF(BFC)  CALL  WRTBFC(  14,  NBFCXE,  XW,  YN,  YS,  ZH,  ZL  , 
SND,NX+1,NY+1,NZ+1,NZ,PRTBFC) 
C     OLD  PRACTICES  RETAINED  FOR  REFERENCE. 
C         IF(SPDATA)  CALL  SPCDAT(IRUN) 
C         IF(BLOCK)  CALL  PORDAT(IRUN) 
IF(GRAPHS)  CALL  SORT(IRUN) 
IF(RESTRT)  GO  TO  902 
DO  901  INDVAR=1,25 

IF(IFIX(FIINIT(INDVAR)+0.1).NE. 10101)  GO  TO  901 
CALL  FLDDAT(IRUN) 
GO  TO  902 

901  CONTINUE 

902  CALL  DATAIO(WRT,10) 

IF(MONITR)  CALL  DATAIOCWRT, -6 ) 
999   CONTINUE 
STOP 
END 
Cxxx    IGEN=1  SO  BFCXYZ  NOT  REQUIRED. 
Cxxx    COMMENT  OUT  BOTH  VERSIONS. 
c 

SUBROUTINE  BFCXYZ  (NXP1 , NYP1 , NZP1 ) 

RETURN 

END 
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CSDIRECTIVEXXMAIN     AMI  LEITNER 

C     LECGRD    LAST  GEO.  NZ=27  NY=18  LAMINAR  FLOW 

r     XFILE  NAME:  MODBFCGD.FTN 

C     XINCLUDE  DED  SUBROUTINES:  THE  MODELS  OF  MAIN,  GROUND  8  STRIDE. 

C     ^DOCUMENTATION:  PHOENICS  INSTRUCTION  MANUAL  (SPRING  1983) 

C     WITH  BODY-FITTED  COORDINATES  INSTRUCTION  SUPPLEMENT 

C      (SUMMER  1984). 

C    ^SATELLITE  FILE  NAME:  MODSTL.FTN 

C0MM0N/ISHIFT/III(57),NFMAX 
C  SET  F-ARRAY  DIMENSION  AS  NEEDED,  8  SET  NFMAX  ACCORDINGLY. 
C      FOR  BFC'S  ALSO  SET  Fl-ARRAY  DIMENSION  AS  NEEDED  ,AND  SET 
C      NF1MAX  ACCORDINGLY. 

COMMON/F0B/FK10000) 

COMMON/NF0B/NF1MAX 

COMMON  F(25000) 

NFMAX=25000 

NF1MAX=10000 

CALL  MAIN1 

STOP 

END 
CSDIRECTIVEXXGROUND 

SUBROUTINE  GROUND( IRN, ICHAP, ISTP, ISWP, IZED, INDVAR) 

INCLUDE  (CMNGUS) 

INCLUDE  (GUSSEQ) 
C       INCLUDE  NMLIST 

LOGICAL  BFC 

EQUIVALENCE  ( LOGIC(20) , BFC) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  1  STARTS: 
c 

C+++++MEANING  OF  SUBROUTINE  ARGUMENTS: 

C      IRN=RUN  NUMBER;  ICHAP=CHAPTER  CALLED;  ISTP=TIME  STEP; 

C      ISWP=SOLUTION  SWEEP;  IZED=Z-SLAB;  INDVAR:  SEE  CHAPTERS  BELOW. 

C+++++USER-INTRODUCED  VARIABLES  &  ARRAYS: 

C      TO  AVOID  CONFLICT  WITH  VARIABLE  NAMES  USED  IN  COMMON,  ALL 

C      VARIABLES  INTRODUCED  BY  THE  USER  SHOULD  HAVE  NAMES  STARTING 

C     WITH  'G'  IF  REAL,  fJ»  IF  INTEGER,  AND  'G'  OR  'J1  IF  LOGICAL. 

C      THUS  GDZ(IZ)  MIGHT  BE  A  Z-INTERVAL  ARRAY; 

C      GW1(IY,IX)  A  2-D  ARRAY  FOR  AXIAL  VELOCITY;  ETC. 

C      USER-GENERATED  SUBROUTINES  SHOULD  BE  NAMED  CORRESPONDINGLY, EG 

C     SUBROUTINE  GVISC(GTEMP,GCNC, GVSC) ,  FOR  COMPUTING  VISCOSITY 

C      FROM  CONCENTRATION  8  TEMPERATURE. 

C+++++GROUND-TO-EARTH  CONNECTING  SUBROUTINES: 


XUSE  GET(NAME,GARRAY,NY,NX)  TO  PUT  VALUES  OF  VARIABLE  NAMED 

•NAME'  INTO  ARRAY  "GARRAY1  DIMENSIONED  GARRAY( NY, NX) . 
XUSE  SET(NAME,IXF,IXL,IYF,IYL,GARRAY,NY,NX)  TO  SET  VARIABLE 
•NAME*  TO  GARRAY(IY,IX)  OVER  THE  REGION:  IXF-IXL  8  IYF-IYL . 
XUSE  PRNSLB(NAME)  TO  PRINT  VARIABLE  'NAME'  OVER  X-Y  PLANE. 
XUSE  ADD(NAME,IXF,IXL,IYF,IYL,TYPE,CM,VM,CVAR,VVAR,NY,NX) 

TO  ADD  SOURCE  TO  VARIABLE  NAMED  'NAME'  (SEE  CHAPTER  5). 
XUSE  READIZ(IZED)  IN  CHAPTERS  1,  2,  8,  8  9  TO  ACCESS  PI,... DM 

8  VOL,...AHDZ.  (SEE  FOOTNOTE  TO  LEGALITY  TABLE) 
XUSE  GET1D(NAME,GARRAY,NDIM)  TO  PUT  VARIABLE  NAMED  'NAME»  IN 
ONE-D  ARRAY  'GARRAY'  DIMENSIONED  NDIM,  THUS: 
CALL  GET1D(NAME,GNX,NX)  FOR  XG,...DXG  8  DIMENSION  GNX(NX); 
CALL  GET1D(NAME,GNY,NY)  FOR  YG,...RV  8  DIMENSION  GNY(NY); 
CALL  GET1D(NAME,GNZ,NZ)  FOR  ZG,...WGRID  8  DIMENSION  GNZ(NZ). 
C+++++LEGALITY  TABLE  FOR  USE  OF  EARTH-CONNECTING  SUBROUTINES: 
C      ENTRIES  IN  TABLE  GIVE  CHAPTERS  IN  WHICH  SUBROUTINES  CAN  BE 
C      USED  FOR  VARIABLES  IN  LEFT-HAND  COLUMN.  (SUBROUTINE 
C      STRIDE  IS  REGARDED  AS  BEING  IN  CHAPTER  3) 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


:  VARIABLE:: 
:  :  : 


GET  8 
PRNSLB 


SET 


ADD 


READIZ  :  GET1D 


:P1  -  RZ  «  : 
:P10  -  RZH: : 
:VOL  -AHDZ: : 
:D1DP  :: 
:D2DP  :: 
:MU1,MU1H  :: 
:EXCO(L,H): : 
.CFP        :: 


3-7 


ALL  : 
10-16: 

ALL  : 

NONE  : 

NONE  : 

5,13-16  : 

NONE  : 

5  : 


>  8 
3 
3 
10 
11 
12 
13 
14 


1,2,8,9:   NONE  : 


NONE 

:    NONE    : 

NONE  : 

NONE 

:  1,2,8,9: 

NONE  : 

NONE 

:    NONE   : 

NONE  : 

NONE 

:    NONE    : 

NONE  : 

NONE 

:    NONE    : 

NONE  : 

NONE 

:    NONE    : 

NONE  : 

NONE 

:    NONE   : 

NONE  : 
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NONE 
NONE 
NONE 


NONE 
NONE 
NONE 


i  NONE  t 
>  NONE  : 
:   ALL   » 


C  iMDT       ::      5     «    15 

C  :HST1,HST2:«   5  8  15   :    16 

C  tXG  -WGRID::     NONE   t   NONE 

c  

C     NOTES  ON  ABOVE  TABLE: 

C     XIN  CHAPTERS  1,  2,  8,  8  9  VARIABLES  PI... DM  8  GEOMETRY 

C     VOL...AHDZ  CAN  BE  ACCESSED  BUT  ONLY  IN  CONJUNCTION  WITH 

C     USE  OF  READIZ,  THUS: 

C     DO  1  IZED=1,NZ 

C     CALL  READIZ(IZED) 

C    1  CALL  GETC...  AS  REQUIRED..) 

C     XGEOMETRY  ACCESSED  BY  READIZ  IS  THAT  AT  INITIAL  TIME. 

C     XD1DP  8  D2DP  ONLY  ACCESSIBLE  IN  UNSTEADY  FLOWS.' 

C+++++GROUND  SERVICE  SUBROUTINES: 

C     XUSE  CONTURCNAME, IPLANE, ILOC, NINT, II, 12, J1,J2, GARRAY, NDIM)  FOR 

C      LINE-PRINTER  PLOTS  OF  CONTOURS.  'NAME'  =  U1,...C4; 

C      «IPLANE'=  XPLANE,  YPLANE,  OR  ZPLANE;  ILOC  SETS  IX,  IY,  OR 

C      IZ  LOCATION  OF  IPLANE;  II,  12,  Jl,  8  J2  SET  FIRST  8  LAST 

C      CELLS  IN  HORIZ.  8  VERT.  ON  PLOT;  GARRAY  IS  1-D  WORKING  ARRAY 

C      OF  DIMENSION  NXXNY,  NXXNZ,  OR  NYXNZ  DICTATED  BY  IPLANE;  8 

C      NDIM  SETS  VALUE  OF  DIMENSION  OF  GARRAY. 

C     XUSE  FLD2DACTITLE, GARRAY, NY, NX)  TO  PRINT  ANY  ARRAY  DIMENSIONED 

C      GARRAY(NY,NX);  SET  'TITLE'  TO  REQUIRED  NAME  (  4  HOLLERITH 

C      CHARACTERS  ONLY). 

C     XUSE  FLD3DA(TITLE, GARRAY, NX, NY, NZ, IPLANE, ILOC)  TO  PRINT  ANY 

C      ARRAY  DIMENSIONED  GARRAYCNX, NY, NZ)  IN  PLANE  SPECIFIED  BY 

C      'IPLANE'  8  'ILOC  AS  FOR  CONTUR  ABOVE;  SET  'TITLE'  AS  FOR 

C      FLD2DA. 

C      VARIABLE  NAMES  FOR  USE  IN  GROUND: 

COMMON/TYPE/CELL, EAST, WEST, NORTH, SOUTH, HIGH, LOW, VOLUME, WALL 

C0MM0N/VAR/P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS, 
SKE,EP,H1,H2,H3,C1,C2,C3,C4,RX,RY,RZ,S1,S2 

C0MM0N/VAR0LD/P10,PP0,U10,U20,V10,V20,W10,W20,R10,R20,RS0, 
&KE0,EP0,H10,H20,H30,C10,C20,C30,C40,RX0,RY0,RZ0,S10,S20 

C0MM0N/VARL0W/P1L,PPL,U1L,U2L,V1L,V2L,W1L,W2L,R1L,R2L,RSL, 
8KEL,EPL,H1L,H2L,H3L,C1L,C2L,C3L,C«L,RXL,RYL,RZL,S1L,S2L 

C0MM0N/VARHI/P1H,PPH,U1H,U2H,V1H,V2H,W1H,W2H,R1H,R2H,RSH, 
SKEH,EPH,HlH,H2H,H3H,ClH,C2H,C3H,C<tH,RXH,RYH,RZH,SlH,S2H 

COMMON/GMTRY/VOL,VOLO,AEAST,ANORTH,AHIGH,AEDX,ANDY,AHDZ 

C0MM0N/PR0P/D1,D2,D1DP,D2DP,MU1,MU1LAM,EXC0,CFP,MDT,HST1,HST2 

COMMON/PRPOLD/DIO, D20 

COMMON/PRPLOW/D1L , D2L , EXCOL 

C0MM0N/PRPHI/D1H,D2H,MU1H,EXC0H 

COMMON/VARNX/XG,XU,DXU,DXG 

COMMON/ VARNY/YG , YV , DYV , DYG , R , RV 

COMMON/VARNZ/ZG, ZW1 , DZW, DZG, WGRID 

COMMON/GDMSCI/XPLANE, YPLANE, ZPLANE, ITNO 

COMMON/GDMSCL/LSLAB,MSLAB,HSLAB,LAMMU 

REAL  NORTH, LOW 

INTEGER     P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS, 
8EP,H1,H2,H3,C1,C2,C3,C4,RX,RY,RZ,S1,S2 

INTEGER     P10,PP0,U10,U20,V10,V20,W10,W20,R10,R20,RS0, 
8EP0,H10,H20,H30,C10,C20,C30,C40,RX0,RY0,RZ0,S10,S20 

INTEGER     P1L,PPL,U1L,U2L,V1L,V2L,W1L,W2L,R1L,R2L,RSL, 
8EPL,H1L,H2L,H3L,C1L,C2L,C3L,C4L,RXL,RYL,RZL,S1L,S2L 

INTEGER     P1H,PPH,U1H,U2H,V1H,V2H,W1H,W2H,R1H,R2H,RSH, 
8EPH,H1H,H2H,H3H,C1H,C2H,C3H,C4H,RXH,RYH,RZH,S1H,S2H 

INTEGER  VOL, VOLO, AEAST, ANORTH, AHIGH, AEDX, ANDY, AHDZ 

INTEGER  D1,D1DP,D2,D2DP,EXC0,CFP,HST1,HST2 

INTEGER  D10,D20,D1L,D2L, EXCOL, D1H, D2H, EXCOH 

INTEGER  XG,XU,DXU,DXG,YG,YV,DYV,DYG,R,RV,ZG,ZW1,DZW, 
8DZG, WGRID 

INTEGER  XPLANE, YPLANE, ZPLANE 

LOGICAL  LSLAB,MSLAB,HSLAB,LAMMU,LSPDA 

EQUIVALENCE  (Ml , Rl ) , (M2, R2) 
C      SATLIT-EQUIVALENT  IRUN: 

EQUIVALENCE  C IRUN, INTGRC 11 ) ) 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  1  ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  1  STARTS: 
C      ARRAYS  (  DIMENSIONED  NY, NX  )  FOR  USE  WITH  'ADD': 

DIMENSION  CVAR(l,l),VVAR(l,l),CM(l,l),VM(l,l),ZERO(l,l) 

DIMENSION  GP(30,1),GH(30,1),GD(30,1),GV(30,1),GW(30,1) 


VAN00730 
VAN00740 
VAN00750 
VAN00760 
VAN00770 
VAN00780 
VAN00790 
VAN00800 
VAN00810 
VAN00820 
VAN00830 
VAN00840 
VAN00850 
VAN00860 
VAN00870 
VAN00880 
VAN00890 
VAN00900 
VAN00910 
VAN00920 
VAN00930 
VAN00940 
VAN00950 
VAN00960 
VAN00970 
VAN00980 
VAN00990 
VAN01O0O 
VAN01010 
VAN01020 
VAN01030 
VAN01040 
VAN01050 
VAN01060 
VAN01070 
VAN01080 
VAN01090 
VAN01100 
VAN01110 
VAN01120 
VAN01130 
VAN011<40 
VAN01150 
VAN01160 
VAN01170 
VAN01180 
VAN01190 
VAN01200 
VAN01210 
VAN01220 
VAN01230 
VAN012<iO 
VAN01250 
VAN01260 
VAN01270 
VAN01280 
VAN01290 
VAN01300 
VAN01310 
VAN01320 
VAN01330 
VAN013^0 
VAN01350 
VAN01360 
VAN01370 
VAN01380 
VAN01390 
VAN01<+00 
VAN01410 
VAN01<i20 
VAN01430 
VANOl'i^O 


59 


FILEs  VANTGRD   FORTRAN   Al 


1  ,GMACH(30,1),GTEMP(30,1),GVISC(30,1),GWH(30,1),GWM(30,1) 

2  ,GKE(30,1),GC3(30,1),GYGC30,1),GXX(30,1),GYY(30,1),GZZ(30,1) 
r     SPECIAL-DATA  ARRAYS  DIMENSIONED  &  DIMENSION  VALUES  SET  HERE: 

DIMENSION  LSPDA(1),ISPDA(1),RSPDA(1) 
C     USER  PLACES  HIS  VARIABLES,  ARRAYS,  EQUIVALENCES  ETC.  HERE. 
EQUIVALENCE  (RAIR, RE(21 ) ) , (GAMA, RE( 22) ) , (GSWP, RE(23)  )  , 
1(GPR,RE(24)),(GTW,RE(25)),(GEMU1,RE(26)),(JEMU1,INTGR(1)) 
DATA  NLSP,NISP,NRSP/1,1,1/ 
DATA  CVAR,VVAR,CM,VM,ZERO/5*0.0/ 
C     USER  PLACES  HIS  DATA  STATEMENTS  HERE. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  1  ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  2  STARTS: 
C     PLEASE  DO  NOT  ALTER,  OR  RE-SET,  ANY  OF  THE  REMAINING 
C      STATEMENTS  OF  THIS  SECTION. 
DATA  NUMCH4  /  0  / 
IF(SPDATA) 
aCALL  RDSPC(IRN,INTGR(12),LSPDA,NLSP,ISPDA,NISP,RSPDA,NRSP) 
CALL  GRDUTY(IRN,ICHAP,IZED,INDVAR) 
IF(BFC)  CALL  BFCGRD( IRN, ICHAP, ISWP, IZED, INDVAR) 
IF(ICHAP.EQ.-5)  GO  TO  10 
IF(ICHAP.LE.0.OR.ICHAP.GT.16)  RETURN 
GO  TO  (100,200,300,4999,500,600,700,800,900,1000,1100,1200, 
SI  300 ,1400, 1500, 16 00), ICHAP 
RETURN 
4999  NUMCH4=  NUMCH4  +  1 

IF  (M0D(NUMCH4,2).EQ.l)  GO  TO  400 
RETURN 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  STANDARD  SECTION  2  ENDS. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  2  STARTS: 

c 

C      CHAPTER  0:  MODIFY  SATLIT  DATA,  AT  START  OF  EACH  IRN. 

c 

10  CONTINUE 
C         IF(.NOT.NAMLST)  RETURN 
C  IF(IRN.EQ.NRUN)  DATFIL= . FALSE. 

C—  READ  SATLIT  DATA  NAMELIST  HERE 

C  CALL  WRIT40(40HENTER  NAMELIST  DATA  FOR  GROUPS 

C  READ(20,G1G24) 

C  CALL  WRIT40C40HENTER  NAMELIST  DATA  FOR  GROUPS 

C  READ(20,G25G42) 
RETURN 


1  TO  24   ) 
25  TO  42  ) 


c 

C      CHAPTER  1:  CALLED  AT  THE  START  OF  EACH  TIME  STEP. 

C      SET  'DT»  HERE  WHEN  TLAST  SET  NEGATIVE  IN  BLOCK  DATA. 

C      'ATIME  +  DT'  GIVES  THE  END  TIME  OF  THE  CURRENT  TIME  STEP. 

C      NOT  ACCESSED  IF  STEADY, OR  PARABOLIC. 

c 

100  CONTINUE 
RETURN 
c 

C     CHAPTER  2:  CALLED  AT  THE  START  OF  EACH  SWEEP. 

c 

200  CONTINUE 

RETURN 
c 

C      CHAPTER  3:  CALLED  AT  THE  START  OF  EACH  SLAB; 

C      NOT  ACCESSED  IF  PARABOLIC,'  BUT  'STRIDE'  IS. 

c 

300  CONTINUE 
RETURN 
c 

C      CHAPTER  4:  CALLED  AT  THE  START  OF  EACH  RE-CALCULATION  OF 

C      VARIABLES  P1,...C4  AT  CURRENT  SLAB.  ITNO=  ITERATION  NUMBER, 
c™ 

400  CONTINUE 
RETURN 
c 

C  CHAPTER  5:  GROUND  CALLED  WHEN  SOURCE  TERM  IS  COMPUTED. 

C  INDVAR  GIVES  DEPENDENT  VARIABLE  IN  QUESTION  IE.  U1,...C4. 

C  TO  ADD  SOURCE  TO  DEPENDENT  VARIABLE  CKSAY)  FOR  IX  =  IXF,IXL 

C  AND  IY=IYF,IYL  INSERT  STATEMENT: 

C  IF(INDVAR.EQ.Cl) 
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VAN01860 
VAN01870 
VAN01880 
VAN01890 
VAN01900 
VAN01910 
VAN01920 
VAN01930 
VAN01940 
VAN01950 
VAN01960 
VAN01970 
VAN01980 
VAN01990 
VAN02000 
VAN02010 
VAN02020 
VAN02030 
VAN02040 
VAN02050 
VAN02060 
VAN02070 
VAN02080 
VAN02090 
VAN02100 
VAN02110 
VAN02120 
VAN02130 
VAN02140 
VAN02150 
VAN02160 
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FILEt  VANTGRD   FORTRAN   Al 

C    8CALL  ADD(INDVAR,IXF,IXL,IYF,IYL,TYPE,CM,VM,CVAR,VVAR,NY,NX)  VAN02170 

C      NOTES  ON  'ADD't  VAN02180 

C    XS0URCE=  (CVAR(IY,IX)+AMAXH0.0,MASFL0))x(VVAR(IY,IX)-PHI),  VAN02190 

C     WHERE  'PHI'  IS  IN-CELL  VALUE  OF  VARIABLE  IN  QUESTION.  VAN02200 

C    x«MASFLO'=  CM(IY,IX)X(VM(IY,IX)-P),  VAN02210 

C     WHERE  *P'  IS  THE  IN-CELL  PRESSURE.  VAN02220 

C     XFOR  INDVAR=  Ml,  OR  =  M2,  SOURCE  ADDED  IS  'MASFLO'  ONLY,  VAN02230 

C      EXCEPT  FOR  ONEPHS=.F.  8  MASFLO  <  0.0  (IE.  OUTFLOW)  WHEN  VAN02240 

C     CM(IY,IX)  IS  MULTIPLIED  BY  R1XD1  (FOR  Ml)  8  R2*D2  (FOR  M2).  VAN02250 

C     XBOTH  'CVAR»  8  'CM'  ARE  MUTLIPLIED  BY  CELL-GEOMETRY  QUANTITY  VAN02260 

C      DICTATED  BY  SETTING  OF  'TYPE'  (=CELL,  EAST  AREA, .. VOLUME) .  VAN02270 

C     XTYPE-SPECIFIED  AREAS  ARE  CALCULATED  AS  IF  BLOCKAGE  ABSENT,  VAN02280 

C      BUT  'VOLUME1  WITH  ACCOUNT  FOR  ITS  PRESENCE.  VAN02290 
C     XFOR  ALL  SOLVED  VARIABLES,  INCLUDE  DING  Ml  (  S  M2  WHEN  ONEPHS=F),   VAN02300 

C      IF  'CM'>  0.0  CALL  'ADD';  FOR  Ml  8  M2  ALTHOUGH  'CVAR'  8  'WAR'  VAN02310 

C      HAVE  NO  SIGNIFICANCE  THEY  MUST  BE  ENTERED  AS  ARGUMENTS.  VAN02320 

C     x'CVAR',  'WAR',  'CM*  8  'VM'  MUST  BE  DIMENSIONED  NY, NX.  VAN02330 

c VAN02340 

500  CONTINUE  VAN02350 

RETURN  VAN02360 

c VAN0237  0 

C      CHAPTER  6:  CALLED  AT  THE  END  OF  EACH  VARIABLE-RECALCULATION  VAN02380 

C      CYCLE  COMMENCED  AT  CHAPTER  4.  ITNO  =  ITERATION  NUMBER.  VAN02390 

C VAN02400 

600  CONTINUE  VAN02410 

RETURN  VAN02420 

c VAN02430 

C      CHAPTER  7:  CALLED  AT  END  OF  EACH  SLAB-WISE  CALCULATION.  VAN02440 

c VAN02450 

700  CONTINUE  '  VAN02460 
IF(FLOATdSWP).LT.GSWP)  RETURN  VAN02470 
CALL  GET(P1,GP,NY,NX)  VAN02480 
CALL  GET(H1,GH,NY,NX)  VAN02490 
CALL  GET(D1,GD,NY,NX)  VAN02500 
CALL  GET(V1,GV,NY,NX)  VAN02510 
CALL  GET(W1,GW,NY,NX)  VAN02520 
CALL  GET(KE,GKE,NY,NX)  VAN02530 

C      CALL  GET1D(YG,GYG,NY)  VAN02540 

CALL  GRED1(39,IZED,GYG,NY,NX)  VAN02550 

CALL  GRED3(57,IZED,GXX,GYY,GZZ,NY,NX)  VAN02560 

GCP=RAIR/(1.-1/GAMA)  VAN02570 

DO  701  1=1, NY  VAN02580 

GSON  =  SQRT(GAMAXGP(I,l)/GD(I,D)  VAN02590 

GAV=SQRT(GV(I,1)XX2+GW(I,1)XK2)  VAN026  00 

GMACH(I,l)=GAV/GSON  VAN02610 

701  GTEMP(I,1)=GP(I,1)/GD(I,1)/RAIR  VAN02620 
C  701  GTEMP(I*l)=(GH(I,l)-GW(I,l)xx2/2.-GV(I,l)x*2/2.)/GCP  VAN02630 

CALL  SET(C1,1,NX,1,NY,GMACH,NY,NX)  VAN02640 

CALL  SET(C2,1,NX,1,NY,GTEMP,NY,NX)  VAN02650 

C CALCULATE  DY1  CF  ST  H(CONVECTIVE  COEF.)  Q  TAU  TR  VAN02660 

IFUEMU1.NE.2)  GOTO  702  VAN02670 

C TURBULENT  VALUES  VAN02680 

GCF=2./GW(NY,1)XX2XGKE(1,1)/3.33XGD(1,1)/GD(NY,1)  VAN02690 

C7     GCF=GCFXGD(NY,1)/GD(1,1)XGTEMP(NY,1)/GTEMP(1,1)XGP(1,1)/GP(NY,1)   VAN027  00 

GST=GCF/2./GPRXM.666  VAN02710 

GHH=GD(NY,1)XGCPXGW(NY,1)XGST  VAN02720 

GR=GPRXX.333  VAN02730 

GTR=GTEMP(NY,1)X(1.+GRX(GAMA-1.)/2.XGMACH(NY,1)XX2)  VAN02740 

C     1(1.+(GAMA-1.)/2.XGMACH(NY,1)X*2)  VAN02750 

GQ=GHHX(GTR-GTW)  VAN02760 

GOTO  703  VAN02770 

C LAMINAR  VALUES  VAN02780 

702  CONTINUE  VAN02790 
IF(JEMUl.EQ.-l)  GEMU1=GVISC(1,1)  VAN02800 
GQ=GEMU1/GPRX(GH(1,1)-GTWXGCP)/GYG(1,1)  VAN02810 
GR=GPRXX.5  VAN02820 
GTR=GTEMP(NY,1)X(1.+GRX(GAMA-1 . )/2 .XGMACH(NY, 1)XX2)  VAN02830 

C     1(1.+(GAMA-1.)/2.XGMACH(NY,1)XX2)  VAN02840 

GHH=GQ/(GTR-GTW)  VAN02850 

GST=GHH/(GD(NY,1)XGW(NY,1)XGCP)  VAN0286  0 

GTAU=GEMU1XGW(1,1)/GYG(1,1)  VAN0287  0 

GCF=GTAUX2./(GD(NY,1)XGW(NY,1)XX2)  VAN0288  0 


61 


FILE:  VANTGRD   FORTRAN   Al 


703   GC3C1,1)=GYG(1,1) 
GC3C2,1)=GCF 
GC3C3,1)=GST 
GC3C4,l)=GCF/2./GST 
GC3C5,1)=GHH 
GC3(6,1)=GQ 
GC3C7,1)=GTAU 
GC3(8,1)=GTR 
GC3C9,1)=GTR-GTW 

GC3C10,l)=GDCNY,l)xGW(NY,l)XGZZCl,l)/GEMUl 
GC3C11,1)=GZZC1,1) 
GC3C12,1)=GEMU1 

GC3C13,l)=GDCNY,l)XGWCNY,l)XGYG(l,l)/GEMUlxSQRTCABS(GCF/2.)) 
CALL  SETCC3,1,NX,1,NY,GC3,NY,NX) 

RETURN 
c — 

C      CHAPTER  8:  CALLED  AT  THE  END  OF  EACH  SWEEP; 

C      NOT  ACCESSED  IF  PARABOLIC. 

c 

800  CONTINUE 

RETURN 
c 

C     CHAPTER  9:  CALLED  AT  THE  END  OF  EACH  TIME  STEP; 

C      NOT  ACCESSED  IF  PARABOLIC. 

c 

900  CONTINUE 

RETURN 
c 

C  CH'^TER  I0«  SET  PHASE  !■  DENSITY  HERE  WHEN  IRH01=-1  IN  DATA. 

C  SET  CURRENT-Z  'SLAB'  DENSITY,  Dl,  IF  MSLAB=.T., 

C  EG.   IFCMSLAB)  CALL  SET( Dl , 1 , NX, 1 , NY, GDI , NY, NX) . 

C  SET  NEXT  LARGER-Z  'SLAB1  DENSITY,  D1H,  IF  HSLAB=.T.  &  PARAB=F 

C  EG.   IF(HSLAB)  CALL  SET( D1H, 1 , NX, 1 , NY,GD1H, NY, NX) . 

C  SET  D(LN(D1))/DP  (IE.  D1DP)  FOR  UNSTEADY  FLOW, 

C      EG.   IF(MSLAB)  CALL  SET( D1DP, 1 , NX, 1 , NY,GD1DP, NY, NX) . 

c 


1000 


101 


102 


CONTINUE 

IF  (MSLAB) 

JP1=P1H 

JH1=H1H 

JD1=D1H 

JW1=W1H 

JV1=V1H 

GO  TO  102 

JP1=P1 

JH1=H1 

JD1=D1 

JW1=W1 

JV1=V1 

CALL  GETCJ 

CALL  GETCJ 

CALL  GETCJ 

CALL  GETCJ 

IFCIZED.EQ 

IFCIZED.EQ 


GO  TO  101 


10^ 

115 
103 


DO  103  IX 

DO  103  IY 

IFCHSLAB) 

GWA=(GW(IY 

GWMCIY,IX) 

GOTO  115 

GWA=CGWCIY 

GWHCIY,IX) 

GHS=GHCIY, 

IFCGHS.LE. 

GDCIY,IX)= 

GOTO  113 


PI, GP, NY, NX) 
HI, GH, NY, NX) 
W1,GW,NY,NX) 
V1,GV,NY,NX) 
.1)  GOTO  105 
.NZ)  GOTO  109 

IZED=2,NZ-1 

I,  NX 

1,NY 

GOTO  104 

,IX)+GWMCIY,IX))/2. 

=GWCIY,IX) 

,IX)+GWHCIY,IX))/2. 
=GWCIY,IX) 

IX)-CGWAX*2+GVCIY,IX)XX2)/2 
1.E5)  GHS=1.E5 
GPC IY, IX)/C 1-1/GAMA)/GHS 


IZED=1 

105  DO  106  IX=1,NX 

DO  106  IY=1,NY 

GHS=GHCIY,IX)-CGWCIY,IX)xx2+GV(IY,IX)xx2)/2 


VAN02890 
VAN02900 
VAN02910 
VAN02920 
VAN02930 
VAN02940 
VAN02950 
VAN02960 
VAN02970 
VAN02980 
VAN02990 
VAN03000 
VAN03010 
VAN03020 
VAN03030 
VAN0304.0 
VAN03050 
VAN03060 
VAN03070 
VAN03080 
VAN03090 
VAN03100 
VAN03110 
VAN03120 
VAN03130 
VAN03140 
VAN03150 
VAN03160 
VAN03170 
VAN03180 
VAN03190 
VAN03200 
VAN03210 
VAN03220 
VAN03230 
VAN03240 
VAN03250 
VAN03260 
VAN03270 
VAN03280 
VAN03290 
VAN03300 
VAN03310 
VAN03320 
VAN03330 
VAN03340 
VAN03350 
VAN03360 
VAN03370 
VAN 03380 
VAN03390 
VAN03400 
VAN03410 
VAN03420 
VAN03430 
VAN03440 
VAN03450 
VAN03460 
VAN03470 
VAN03480 
VAN03490 
VAN03500 
VAN03510 
VAN03520 
VAN03530 
VAN03540 
VAN03550 
VAN03560 
VAN03570 
VAN0358<r 
VAN03590 
VAN03600 
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107 
106 


109 


111 


112 

110 


IF(GHS.LE.1.E5)  GHS=1.E5 
GD(IY,IX)=  GP(IY,IX)/(1-1/GAMA)/GHS 
IF(HSLAB)  GOTO  107 
GWM(IY,IX)=GW(IY,IX) 
GOTO  106 

GWH(IY,IX)=GW(IY,IX) 
CONTINUE 
GOTO  113 
IZED=NZ 

DO  110  IX=1,NX 

DO  110  IY=1,NY 

IF(HSLAB)  GOTO  111 

GHS=GH(IY,IX)-(GWM(IY,IX)XX2+GV(IY,IX)XX2)/2. 

IF(GHS.LE.1.E5)  GHS=1.E5 

GWM(IY,IX)=GW(IY,IX) 

GOTO  112 

GHS=GH(IY,IX)-(GWH(IY,IX)X*2+GV(IY,IX)xx2)/2. 

IF(GHS.LE.1.E5)  GHS=1.E5 

GWH(IY,IX)=GW(IY,IX) 

GD(IY,IX)=  GP(IY,IX)/(1-1/GAMA)/GHS 

CONTINUE 


C- 
C 

C 
C 

c 
c 
c 
c 

c- 


113  CONTINUE 

CALL  SET(JD1,1,NX,1,NY,GD,NY,NX) 
RETURN 


C- 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C- 


CHAPTER  11:  SET  PHASE  2  DENSITY  HERE  WHEN  IRH02=-1  IN  DATA. 

SET  CURRENT-Z  'SLAB'  DENSITY,  D2,  IF  MSLAB=.T., 

EG.   IF(MSLAB)  CALL  SET( D2, 1 , NX, 1 , NY,GD2, NY, NX) . 

SET  NEXT  LARGER-Z  'SLAB'  DENSITY,  D2H,  IF  HSLAB=.T.  &  PARAB=F 

EG.   IF(HSLAB)  CALL  SET( D2H, 1 , NX, 1 , NY,GD2H, NY, NX) . 

SET  D(LN(D2))/DP  FOR  UNSTEADY  FLOW, 

EG.   IF(MSLAB)  CALL  SET( D2DP, 1 , NX, 1 , NY, GD2DP, NY, NX) . 

1100  CONTINUE 
RETURN 

CHAPTER  12s  SET  PHASE  1  VISCOSITY  HERE  WHEN  IEMU1=-1  IN  DATA. 

SET  CURRENT-Z  'SLAB'  VISCOSITY  (MUD,  IF  MSLAB=.T., 

EG.   IF(MSLAB)  CALL  SET(MU1 , 1 , NX, 1 , NY, GVISC, NY, NX) . 

SET  NEXT  LARGER-Z  'SLAB'  VISC.  (MU1H),  IF  HSLAB=.T.  &  PARAB=F 

EG.   IF(HSLAB)  CALL  SETCMU1H, 1 , NX, 1 , NY,GVSCH, NY, NX) . 

CHAPTER  ALSO  ACCESSED  WHEN  EMULAM=-1.0  IN  DATA,  SO  THAT  THE 
LAMINAR  VISCOSITY  WHICH  APPEARS  IN  WALL  FUNCTIONS  8.  IN  THE 
KE-EP  TURBULENCE  MODEL  (IEMU1=2)  MAY  BE  SET  NON-CONSTANT. 
SET  CURRENT-Z  'SLAB'  VALUE  (MU1LAM)  WHEN  LAMMU=.T., 
EG.   IF(LAMMU)  CALL  SET(MU1LAM, 1 , NX, 1 , NY,GVSCL , NY, NX) . 


1200 


122 
123 


121 
C121 


CONTINUE 

GCP=RAIR/(1.-1/GAMA) 

IF  (HSLAB)   GOTO  122 

CALL  GET(H1,GH,NY,NX) 

CALL  GET(W1,GW,NY,NX) 

CALL  GET(V1,GV,NY,NX) 

GOTO  123 

CALL  GET(H1H,GH,NY,NX) 

CALL  GET(W1H,GW,NY,NX) 

CALL  GETCV1H,GV,NY,NX) 

CONTINUE 

DO  121  IX=1,NX 

DO  121  IY=1,NY 

GTMP=(GH(IY,IX)-GW(IY,IX)xx2/2.-GV(IY,IX)xx2/2.)/GCP 

IFCGTMP.LT. 150.)   GTMP=150. 

GVISCCIY,IX)=1.716E-05X(GTMP/273.)XX0.666 

IF(GVISC(IY,IX).LE. .8E-5)   GVISCC IY, IX) = . 8E-5 

IF  (MSLAB)  CALL  SET(MU1 , 1 , NX, 1 , NY, GVISC, NY, NX) 

IF  (HSLAB)  CALL  SET(MU1H, 1 , NX, 1 , NY, GVISC, NY, NX) 

IF  (LAMMU)  CALL  SET(MU1LAM, 1 , NX, 1 , NY, GVISC, NY, NX) 

RETURN 

CHAPTER  13:  SET  EXCHANGE  COEFFICIENT  (E.C.)  FOR  VARIABLE 


VAN03610 
VAN03620 
VAN03630 
VAN03640 
VAN03650 
VAN03660 
VAN03670 
VAN03680 
VAN03690 
VAN03700 
VAN03710 
VAN03720 
VAN03730 
VAN03740 
VAN03750 
VAN03760 
VAN03770 
VAN03780 
VAN03790 
VAN03800 
VAN03810 
VAN03820 
VAN03830 
VAN03840 
VAN03850 
VAN03860 
VAN03870 
VAN03880 
VAN03890 
VAN03900 
VAN03910 
VAN03920 
VAN03930 
VAN03940 
VAN03950 
VAN03960 
VAN03970 
VAN03980 
VAN03990 
VAN04000 
VAN04010 
VAN04020 
VAN04030 
VAN04040 
VAN04050 
VAN04060 
VAN04070 
VAN04080 
VAN04090 
VAN04100 
VAN04110 
VAN04120 
VAN04130 
VAN04140 
VAN04150 
VAN04160 
VAN04170 
VAN04180 
VAN04190 
VAN04200 
VAN04210 
VAN04220 
VAN04230 
VAN04240 
VAN04250 
VAN04260 
VAN04270 
VAN04280 
VAN04290 
VAN04300 
VAN04310 
VAN04320 
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C     INDVAR  WHEN  SIGMAC INDVAR)=-1 . 0  IN  DATA.  VAN04330 

C     SET  CURRENT-Z  'SLAB'  E.C.  (EXCO)  IF  MSLAB=.T.,  VAN04340 

C     EG.   IF(MSLAB)  CALL  SET( EXCO, 1 , NX, 1 , NY,GEXCO, NY, NX) .  VAN04350 

C     SET  NEXT  SMALLER-Z  'SLAB'  E.C.  (EXCOL)  IF  LSLAB=.T.,  VAN04360 

C     EG.   IF(LSLAB)  CALL  SET( EXCOL , 1 , NX, 1 , NY, GEXCOL , NY, NX) .  VAN04370 

C     SET  NEXT  LARGER-Z   "SLAB1  E.C.  (EXCOH)  IF  HSLAB=.T.,  VAN04380 

C      EG.   IF(HSLAB)  CALL  SET( EXCOH, 1 , NX, 1 , NY, GEXCOH, NY, NX) .  VAN04390 

C     NOTE:  FOR  MSLAB,  INDVAR=U1 , . .C4;  FOR  LSLAB,  INDVAR=U1L , . . C4L  VAN04400 

C     &  FOR  HSLAB,  INDVAR=U1H, . .C4H .  IF  PARAB=.T.  SET  MSLAB  ONLY.  VAN04410 

c VAN04420 

1300  CONTINUE  VAN04430 

RETURN  VAN04440 

c VAN04450 

C     CHAPTER  14:  SET  INTER-PHASE  FRICTION  COEFFICIENT  (CFP)  HERE  VAN04460 

C     WHEN  ICFIP  =  -1  IN  DATA;  ITS  UNITS  =  FORCE  /  (CELL  x  RELATIVE  VAN04470 

C      SPEED  OF  PHASES).  VAN04480 

c VAN0449  0 

1400  CONTINUE  VAN04500 

RETURN  VAN04510 

c VAN04520 

C      CHAPTER  15:  SET  INTER-PHASE  MASS-TRANSFER  RATE  PER  CELL  (MDT)  VAN04530 

C      HERE  WHEN  IMDOT  =  -1  IN  DATA.  VAN04540 

c VAN04550 

1500  CONTINUE  VAN04560 

RETURN  VAN04570 

c VAN04580 

C      CHAPTER  16:  SET  HERE  PHASE  1  &  2  SATURATION  ENTHALPIES  VAN04590 

C      (  HST1  &  HST2)  WHEN  IHSAT  =  -1  IN  DATA.  VAN04600 

c VAN04610 

1600  CONTINUE  VAN04620 

RETURN  VAN04630 

END  VAN04640 
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